Method for simulating chatter-free milled surface topography

ABSTRACT

Method for simulating chatter-free milled surface topography includes: performing dynamics modeling on milling system and building delay differential equation with multiple delays to represent dynamic model; constructing state transition matrix between two adjacent cutting tool rotation periods by extending GRK method; predicting stable milling parameter zones based on Floquet theory; calculating vibration displacements on discretized time nodes during one cutting tool rotation period by means of fixed mapping point theorem; constructing motion trajectories of cutting tool edges in normal and feed directions with regard to milled surface of workpiece; predicting surface topography of workpiece by using spline interpolation to densify motion trajectories of cutting tool edges that participate in generating final surface of workpiece; calculating surface location error of milled surface based on motion trajectories of cutting tool edges in normal direction with regard to milled surface of workpiece; calculating surface roughness based on the simulated surface topography of workpiece.

FIELD OF THE INVENTION

The invention relates to mechanical machining, and more specifically, to a method for simulating chatter-free milled surface topography.

BACKGROUND

Topography of milled surfaces is generated by complicated interaction between cutting tool and workpiece during machining process. In order to achieve qualified milled surface, it is a prerequisite to avoid regenerative milling chatter, which is a kind of self-excited strong vibrations. Nevertheless, chatter-free milling does not mean high-quality milling. Surface location error and surface roughness are closely related to tool geometry (pitch angles, helix angles, number of teeth etc.), cutting conditions (tool runout) and process parameters (spindle rotation speed, feed per tooth, radial/axial depth of cut etc.).

Currently, initial value based time-domain simulation method is the only way to simultaneously predicting chatter stability, surface location error and surface roughness for milling operations, as reported in REF. 1 (Schmitz, T. L., and Smith, K. S., 2009, Machining Dynamics, Springer US, Boston, Mass.). However, the initial value based time-domain simulation method suffers from low computational efficiency, which prevents its wide application in industry.

REF 2 (Niu, J., Ding, Y, Zhu, L., and Ding, H., 2014, Runge-Kutta Methods for a Semi Analytical Prediction ofMilling Stability, Nonlinear Dynamics, 76(1), pp. 289-304) presented a semi-analytical method (i.e., Generalized Runge-Kutta (GRK) method) for stability analysis of milling operations using tools with constant pitch and constant helix angles, and omitting influences of tool runout, which is of high computational accuracy and efficiency. REF. 3 (Niu, J., Ding, Y, Zhu, L., and Ding, H., 2017, Mechanics and Multi-Regenerative Stability of Variable Pitch and Variable Helix Milling Tools Considering Runout, Int. J. Mach. Tools Manuf., 123, pp. 129-145) extended the GRK method in REF. 2 to predicting stability of milling operations using variable pitch and variable helix tools and considering effect of runout. Nevertheless, both REF. 2 and REF. 3 omitted steady cutting force term that does not affect chip regeneration. As a result, the methods in REF. 2 and REF. 3 cannot be used directly for predicting surface location error and surface roughness of milled surface, and cannot be used for simulating surface topography either.

Chinese patents with publication numbers CN102490081B, CN103713576B and CN108515217B also proposed different methods to simulate milled surface topography. However, only kinematic motions were taken into account in these patents, and influences of cutting vibrations were neglected.

In conclusion, it is necessary to propose a method for simulating chatter-free milled surface topography, which can predict chatter stability, surface location error and surface roughness simultaneously. It is of great significance for avoidance of cutting chatter and improvement of cutting quality.

SUMMARY

This invention discloses a method for simulating chatter-free milled surface topography. Its basic idea is: achieving chatter-free cutting parameters by means of dynamics modeling and stability analysis of milling system; calculating relative vibration displacements between cutting tool and workpiece by extending the GRK method; obtaining motion trajectories of cutting tool edges by means of kinematic synthesis of feed movements along tool path, rotations around spindle axis and dynamic vibrations of milling system; and finally obtaining milled surface topography by means of Boolean operation on envelope of cutting tool edge trajectories and workpiece block.

This invention is implemented by following technical solutions.

A method for simulating chatter-free milled surface topography, comprising steps of:

step 1: performing dynamics modeling on milling system and building a delay differential equation with multiple delays to represent dynamic model;

step 2: constructing state transition matrix between two adjacent cutting tool rotation periods by extending Generalized Runge-Kutta method;

step 3: predicting stable milling parameter zones based on Floquet theory;

step 4: calculating vibration displacements on discretized time nodes during one cutting tool rotation period by means of fixed mapping point theorem;

step 5: constructing motion trajectories of cutting tool edges in normal and feed directions with regard to milled surface of workpiece;

step 6: predicting surface topography of workpiece by using spline interpolation to densify motion trajectories of cutting tool edges that participate in generating final surface of workpiece;

step 7: calculating surface location error of milled surface based on motion trajectories of cutting tool edges in normal direction with regard to milled surface of workpiece;

step 8: calculating surface roughness based on the simulated surface topography of workpiece.

The step 1 comprises sub-steps of:

step 1.1: considering flexibility of both cutting tool and workpiece, variance of pitch angles and helix angles of cutting tool, and teeth runout, perform dynamics modeling on milling system and build a delay differential equation with multiple delays to represent the dynamic model in physical space:

${{M{\overset{¨}{q}(t)}} + {C{\overset{˙}{q}(t)}} + {K{q(t)}}} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{{K_{c}\left( {t,j,k} \right)}\left\lbrack {{q(t)} - {q\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack} + {F_{0}\left( {t,j,k} \right)}} \right\}}}$

where M, C and K are mass, damping and stiffness matrices, respectively; {umlaut over (q)}(t), {dot over (q)}(t) and q(t) are acceleration, velocity and displacement vectors at time t, respectively; K_(c)(t,j,k) is directional coefficient matrix for Tooth k at j-th axial slice at time t; F₀(t,j,k) is cutting force vector for Tooth k at j-th axial slice at time t, which is independent of chip regeneration; Σ_(p=k) _(v) ^(k) τ(j,p) is time delay with regard to cutting element (j,k), Σ is mathematical summation operator, p is a variable for summation, and k_(v) is tooth serial number to start delay calculation; N is number of teeth; and N_(a) is number of axial discretization of tool;

step 1.2: perform modal coordinate transformation on the dynamic model in step 1.1, which can be transformed from physical space into modal space; the modal coordinate transformation equation is:

q(t)=PΓ(t)

where P is the mode shape matrix, and Γ(t) is the modal displacement at time t;

dynamics model in modal space which is represented with delay differential equation with multiple delays turns out to be:

${{M_{\Gamma}{\overset{¨}{\Gamma}(t)}} + {C_{\Gamma}{\overset{.}{\Gamma}(t)}} + {K_{\Gamma}{\Gamma(t)}}} = {\underset{k = 1}{\sum\limits^{N}}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{P^{T}{K_{c}\left( {t,j,k} \right)}{P\left\lbrack {{\Gamma(t)} - {\Gamma\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack}} + {P^{T}{F_{0}\left( {t,j,k} \right)}}} \right\}}}$

where M_(Γ), C_(Γ) and K_(Γ) are modal mass, modal damping and modal stiffness matrices, respectively; {circumflex over (Γ)}(t), {dot over (Γ)}(t) and Γ(t) are modal acceleration, modal velocity and modal displacement vectors, respectively;

for milling system with flexible cutting tool and flexible workpiece, the dynamic model becomes:

${{\begin{bmatrix} M_{\Gamma;T} & \\  & M_{\Gamma;W} \end{bmatrix}\begin{Bmatrix} {{\overset{¨}{\Gamma}}_{T}(t)} \\ {{\overset{¨}{\Gamma}}_{W}(t)} \end{Bmatrix}} + {\begin{bmatrix} C_{\Gamma;T} & \\  & C_{\Gamma;W} \end{bmatrix}\begin{Bmatrix} {{\overset{˙}{\Gamma}}_{T}(t)} \\ {{\overset{˙}{\Gamma}}_{W}(t)} \end{Bmatrix}} + \text{ }{\begin{bmatrix} K_{\Gamma;T} & \\  & K_{\Gamma;W} \end{bmatrix}\begin{Bmatrix} {\Gamma_{T}(t)} \\ {\Gamma_{W}(t)} \end{Bmatrix}}} = \text{ }{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{\begin{bmatrix} P_{T}^{T} \\ {- P_{W}^{T}} \end{bmatrix}{{K_{c}\left( {t,j,k} \right)}\begin{bmatrix} P_{T} & {- P_{W}} \end{bmatrix}}\begin{Bmatrix} {{\Gamma_{T}(t)} - {\Gamma_{T}\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \\ {{\Gamma_{W}(t)} - {\Gamma_{W}\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \end{Bmatrix}} + \text{ }{\begin{bmatrix} P_{T}^{T} \\ {- P_{W}^{T}} \end{bmatrix}{F_{0}\left( {t,j,k} \right)}}} \right\}}}$

where m_(Γ;T), C_(Γ;T) and K_(Γ;T) are modal mass, modal damping and modal stiffness matrices of cutting tool; M_(Γ;W), C_(Γ;W) and K_(Γ;W) are modal mass, modal damping and modal stiffness matrices of workpiece; {umlaut over (Γ)}_(T)(t), {dot over (Γ)}_(T) (t) and Γ_(T)(t) are modal acceleration, modal velocity and modal displacement vectors of cutting tool; {umlaut over (Γ)}_(W)(t), {dot over (Γ)}_(W)(t) and Γ_(W)(t) are modal acceleration, modal velocity and modal displacement vectors of workpiece; P_(T) and P_(W) are mode shape matrices of cutting tool and workpiece, respectively; P_(T) ^(T) and P_(W) ^(T) are transposed matrices of P_(T) and P_(W), respectively;

for milling system with flexible cutting tool and rigid workpiece, the dynamic model becomes:

${{M_{\Gamma;T}{{\overset{¨}{\Gamma}}_{T}(t)}} + {C_{\Gamma;T}{{\overset{˙}{\Gamma}}_{T}(t)}} + {K_{\Gamma;T}{\Gamma_{T}(t)}}} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{P_{T}^{T}{K_{c}\left( {t,j,k} \right)}\left\{ {P_{T}\left\lbrack {{\Gamma_{T}(t)} - {\Gamma_{T}\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack} \right\}} + {P_{T}^{T}{F_{0}\left( {t,j,k} \right)}}} \right\}}}$

for milling system with rigid cutting tool and flexible workpiece, the dynamic model becomes:

${{M_{\Gamma;W}{{\overset{¨}{\Gamma}}_{W}(t)}} + {C_{\Gamma;W}{{\overset{˙}{\Gamma}}_{W}(t)}} + {K_{\Gamma;W}{\Gamma_{W}(t)}}} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}{\left\{ {{P_{W}^{T}{K_{c}\left( {t,j,k} \right)}\left\{ {P_{W}\left\lbrack {{\Gamma_{W}(t)} - {\Gamma_{W}\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack} \right\}} - {P_{W}^{T}{F_{0}\left( {t,j,k} \right)}}} \right\}.}}}$

The step 2 comprises sub-steps of:

step 2.1: perform state space transformation on the dynamic model in modal space by introducing state variable vector x(t)=[(Γ(t))^(T),({dot over (Γ)}(t))^(T)]^(T), and delay differential equation with multiple delays in state space becomes:

${\overset{.}{x}(t)} = {{A{x(t)}} + {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{{B\left( {t,j,k} \right)}\left\lbrack {{x(t)} - {x\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack} + {D\left( {t,j,k} \right)}} \right\}}}}$

where x(t) denotes state variable vector, {dot over (x)}(t) denotes first derivative of x(t) with regard to time t,

${A = \begin{bmatrix} 0 & I \\ {{- M_{\Gamma}^{- 1}}K_{\Gamma}} & {{- M_{\Gamma}^{- 1}}C_{\Gamma}} \end{bmatrix}}\ ,$ ${{B\left( {t,j,k} \right)} = \begin{bmatrix} 0 & 0 \\ {M_{\Gamma}^{- 1}P^{T}{K_{c}\left( {t,j,k} \right)}P} & 0 \end{bmatrix}},$ ${{D\left( {t,j,k} \right)} = \begin{bmatrix} 0 \\ {M_{\Gamma}^{- 1}P^{T}{F_{0}\left( {t,j,k} \right)}} \end{bmatrix}},$

and I is unit matrix;

step 2.2: divide one rotation period of cutting tool T into 2i_(m)+2 equal intervals with a series of discrete nodes {t₀,t₁, . . . ,t_(2i) _(m) ₊₂}; analytical expression of dynamic responses at interval [t_(i),t] with regard to the delay differential equation with multiple delays in step 2.1 is:

${x(t)} = {{e^{A({t - t_{i}})}{x\left( t_{i} \right)}} + \text{ }{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}{\int_{t_{i}}^{t}{\left\{ {{e^{A({t - \xi})}{{B\left( {\xi,j,k} \right)}\left\lbrack {{x(\xi)} - {x\left( {\xi - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack}} + \text{ }{e^{A({t - \xi})}{D\left( {\xi,j,k} \right)}}} \right\} d\xi}}}}}$

where ξ is integration variable.

for brevity, define x_(i):=x(t_(i)), E(t,t_(i)):=e^(A(t−t) i ⁾x(t_(i)), H_(j,k)(t,ξ,x(ξ):=e^(A(t−ξ))B(ξ,j,k)[x(ξ)−x(ξ−τ(j,k,p))] and J_(j,k)(t,ξ):=e^(A(t−ξ))D(ξ,j,k); on sub-interval [t_(2i),t_(2i+2)], x(t) is approximated using Simpson's rule:

$x_{{2i} + 2} = {{E\left( {t_{{2i} + 2},t_{2i}} \right)} + {\frac{h}{3}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{H_{j,k}\left( {t_{{2i} + 2},t_{2i},x_{2i}} \right)} + {4{H_{j,k}\left( {t_{{2i} + 2},t_{{2i} + 1},x_{{2i} + 1}} \right)}} + \text{ }{H_{j,k}\left( {t_{{2i} + 2},t_{{2i} + 2},x_{{2i} + 2}} \right)} + {J_{j,k}\left( {t_{{2i} + 2},t_{2i}} \right)} + {4{J_{j,k}\left( {t_{{2i} + 2},t_{{2i} + 1}} \right)}} + {J_{j,k}\left( {t_{{2i} + 2},t_{{2i} + 2}} \right)}} \right\}}}}}$

on sub-interval [t_(2i),t_(2i+1)], x(t) is approximated using fourth-order Runge-Kutta algorithm:

$x_{{2i} + 1} = {{E\left( {t_{{2i} + 1},t_{2i}} \right)} + {\frac{h}{6}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{H_{j,k}\left( {t_{{2i} + 1},t_{2i},x_{2i}} \right)} + {2{H_{j,k}\left( {t_{{2i} + 1},t_{{2i} + {1/2}},x_{{2i} + {1/2}}} \right)}} + \text{ }{2{H_{j,k}\left( {t_{{2i} + 1},t_{{2i} + {1/2}},x_{{2i} + {1/2}}} \right)}} + {H_{j,k}\left( {t_{{2i} + 1},t_{{2i} + 1},x_{{2i} + 1}} \right)} + {J_{j,k}\left( {t_{{2i} + 1},t_{2i}} \right)} + {2{J_{j,k}\left( {t_{{2i} + 1},t_{{2i} + {1\text{/2}}}} \right)}} + {2{J_{j,k}\left( {t_{{2i} + 1},t_{{2i} + {1/2}}} \right)}} + {J_{j,k}\left( {t_{{2i} + 1},t_{{2i} + 1}} \right)}} \right\}}}}}$

where state variable vector at middle points t_(2i+1/2) is approximated with Lagrange formula:

x _(2i+)½=⅜x _(2i)+¾x _(2i+1)⅛x _(2i+2)

the analytical expression of dynamic responses on sub-interval [t_(2i),t_(2i+1)] is re-organized into:

${{F_{2i}x_{2i}} + {F_{{2i} + 1}x_{{2i} + 1}} + {F_{{2i} + 2}x_{{2i} + 2}}} = {{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left( {{F_{{2i} - {\tau({j,k})}}x_{{2i} - {\tau({j,k})}}} + {F_{{2i} + 1 - {\tau({j,k})}}x_{{2i} + {\tau({j,k})}}} + {F_{{2i} + 2 - {\tau({j,k})}}x_{{2i} + 2 - {\tau({j,k})}}}} \right)}} + {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left( {S_{2i} + S_{{2i} + {1/2}} + S_{{2i} + 1}} \right)}}}$ where ${F_{2i} = {{- e^{Ah}} - {\frac{h}{6}e^{Ah}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{2i},j,k}}}} - {{\frac{3}{8} \cdot \frac{4h}{6}}e^{{Ah}/2}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + {1/2}},j,k}}}}}},$ ${F_{{2i} + 1} = {I - {{\frac{3}{4} \cdot \frac{4h}{6}}e^{{Ah}/2}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + {1/2}},j,k}}}} - {\frac{h}{6}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + 1},j,k}}}}}},$ ${F_{{2i} + 2} = {{\frac{1}{8} \cdot \frac{4h}{6}}e^{{Ah}/2}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + {1/2}},j,k}}}}},$ ${F_{{2i} - {\tau({j,k,p})}} = {{{- \frac{h}{6}}e^{Ah}B_{{2i},j,k}} - {{\frac{3}{8} \cdot \frac{4h}{6}}e^{{Ah}/2}B_{{{2i} + {1/2}},j,k}}}},$ ${F_{{2i} + 1 - {\tau({j,k,p})}} = {{{{- \frac{3}{4}} \cdot \frac{4h}{6}}e^{{Ah}/2}B_{{{2i} + {1/2}},j,k}} - {\frac{h}{6}B_{{{2i} + 1},j,k}}}},$ ${F_{{2i} + 2 - {\tau({j,k,p})}} = {{\frac{1}{8} \cdot \frac{4h}{6}}e^{{Ah}/2}B_{{{2i} + {1/2}},j,k}}},$ ${S_{2i} = {\frac{h}{6}e^{A \cdot h}D_{{2i},j,k}}},$ ${S_{{2i} + {1/2}} = {{\frac{h}{6} \cdot 4}e^{{A \cdot h}/2}D_{{{2i} + {1/2}},j,k}}},$ ${S_{{2i} + 1} = {\frac{h}{6}D_{{{2i} + 1},j,k}}},$

and h is discretization step;

the analytical expression of dynamic responses on sub-interval [t_(2i),t_(2i+2)] is re-organized into:

${{G_{2i}x_{2i}} + {G_{{2i} + 1}x_{{2i} + 1}} + {G_{{2i} + 2}x_{{2i} + 2}}} = {{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left( {{G_{{2i} - {\tau({j,k})}}x_{{2i} - {\tau({j,k})}}} + {G_{{2i} + 1 - {\tau({j,k})}}x_{{2i} + 1 - {\tau({j,k})}}} + {G_{{2i} + 2 - {\tau({j,k})}}x_{{2i} + 2 - {\tau({j,k})}}}} \right)}} + {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left( {T_{2i} + T_{{2i} + 1} + T_{{2i} + 2}} \right)}}}$ where ${G_{2i} = {{- e^{2Ah}} - {\frac{h}{3}e^{2Ah}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{2i},j,k}}}}}},$ ${G_{{2i} + 1} = {{- \frac{4h}{3}}e^{Ah}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + 1},j,k}}}}},$ ${G_{{2i} + 2} = {I - {\frac{h}{3}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + 2},j,k}}}}}},$ ${G_{{2i} - {\tau({j,k,p})}} = {{- \frac{h}{3}}e^{2Ah}B_{{2i},j,k}}},$ ${G_{{2i} + 1 - {\tau({j,k,p})}} = {{- \frac{4h}{3}}e^{Ah}B_{{{2i} + 1},j,k}}},$ ${G_{{2i} + 2 - {\tau({j,k,p})}} = {{- \frac{h}{3}}B_{{{2i} + 2},j,k}}},$ ${T_{2i} = {\frac{h}{3}e^{{A \cdot 2}h}D_{{2i},j,k}}},$ $T_{{2i} + 1} = {{\frac{h}{3} \cdot 4}e^{A \cdot h}D_{{{2i} + 1},j,k}}$ and ${T_{{2i} + 2} = {\frac{h}{3}D_{{{2i} + 2},j,k}}};$

step 2.3: based on the formulae in step 2.2, discrete mapping relationship between two adjacent rotation periods of cutting tool [−T,0] and [0,T] is established:

P₁y_([0, T]) = Qy_([−T, 0]) + P₂y_([0, T]) + z_([0, T]) where ${y_{\lbrack{0,T}\rbrack} = \begin{bmatrix} x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\  \vdots \\ x_{{2i_{m}} + 1} \\ x_{{2i_{m}} + 2} \end{bmatrix}},$ ${y_{\lbrack{{- T},0}\rbrack} = \begin{bmatrix} x_{0 - T} \\ x_{1 - T} \\ x_{2 - T} \\ x_{3 - T} \\ x_{4 - T} \\  \vdots \\ x_{{2i_{m}} + 1 - T} \\ x_{{2i_{m}} + 2 - T} \end{bmatrix}},$ ${z_{\lbrack{0,T}\rbrack} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\begin{bmatrix} 0 \\ {S_{0} + S_{1/2} + S_{1}} \\ {T_{0} + T_{1} + T_{2}} \\ {S_{2} + S_{2 + {1/2}} + S_{3}} \\ {T_{2} + T_{3} + T_{4}} \\  \vdots \\ {S_{2i_{m}} + S_{{2i_{m}} + {1/2}} + S_{{2i_{m}} + 1}} \\ {T_{2i_{m}} + T_{{2i_{m}} + 1} + T_{{2i_{m}} + 2}} \end{bmatrix}}}},$ ${P_{1} = \begin{bmatrix} I & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ F_{0} & F_{1} & F_{2} & 0 & 0 & \cdots & 0 & 0 & 0 \\ G_{0} & G_{1} & G_{2} & 0 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & F_{2} & F_{3} & F_{4} & \cdots & 0 & 0 & 0 \\ 0 & 0 & G_{2} & G_{3} & G_{4} & \cdots & 0 & 0 & 0 \\  \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & 0 & \cdots & F_{2i_{m}} & F_{{2i_{m}} + 1} & F_{{2i_{m}} + 2} \\ 0 & 0 & 0 & 0 & 0 & \cdots & G_{2i_{m}} & G_{{2i_{m}} + 1} & G_{{2i_{m}} + 2} \end{bmatrix}},$ $Q = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\begin{bmatrix} 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 & {I/{NN}_{a}} \\ 0 & \cdots & F_{- {\tau({j,k,p})}} & F_{1 - {\tau({j,k,p})}} & F_{2 - {\tau({j,k,p})}} & \cdots & 0 & 0 & 0 \\ 0 & \cdots & G_{- {\tau({j,k,p})}} & G_{1 - {\tau({j,k,p})}} & G_{2 - {\tau({j,k,p})}} & \ddots & \vdots & 0 & 0 \\  \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & 0 & 0 & \cdots & F_{{2i} - {\tau({j,k,p})}} & F_{{2i} + 1 - {\tau({j,k,p})}} & F_{{2i} + 2 - {\tau({j,k,p})}} \\ 0 & \cdots & 0 & 0 & 0 & \cdots & G_{{2i} - {\tau({j,k,p})}} & G_{{2i} + 1 - {\tau({j,k,p})}} & G_{{2i} + 2 - {\tau({j,k,p})}} \\  \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \end{bmatrix}}}$ and $P_{2} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}{\begin{bmatrix} 0 & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 \\  \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \cdots & 0 \\ F_{{2i} + 2 - {\tau({j,k,p})}} & F_{{2i} + 3 - {\tau({j,k,p})}} & F_{{2i} + 4 - {\tau({j,k,p})}} & \cdots & 0 & 0 & 0 & \cdots & 0 \\ G_{{2i} + 2 - {\tau({j,k,p})}} & G_{{2i} + 3 - {\tau({j,k,p})}} & G_{{2i} + 4 - {\tau({j,k,p})}} & \ddots & 0 & 0 & 0 & \cdots & 0 \\ 0 & 0 & \vdots & \ddots & 0 & 0 & 0 & \cdots & 0 \\  \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & F_{{2i_{m}} - {\tau({j,k,p})}} & F_{{2i_{m}} + 1 - {\tau({j,k,p})}} & F_{{2i_{m}} + 2 - {\tau({j,k,p})}} & \cdots & 0 \\ 0 & 0 & 0 & \cdots & G_{{2i_{m}} - {\tau({j,k,p})}} & G_{{2i_{m}} + 1 - {\tau({j,k,p})}} & G_{{2i_{m}} + 2 - {\tau({j,k,p})}} & \cdots & 0 \end{bmatrix}.}}}$

The step 3 is:

construct transition matrix Φ of the milling system:

$\Phi = \left\{ \begin{matrix} {{\left( {P_{1} - P_{2}} \right)^{- 1}Q},} & {{{if}{❘{P_{1} - P_{2}}❘}} \neq 0} \\ {{\left( {P_{1} - P_{2}} \right)^{\dagger}Q},} & {{{if}{❘{P_{1} - P_{2}}❘}} = 0} \end{matrix} \right.$

where |P₁−P₂| denotes determinant of matrix (P₁−P₂); † denotes Moore-Penrose pseudoinverse of matrix (P₁−P₂); based on Floquet theory, stability of milling system depends on eigenvalue of Φ: it is stable if all modules of eigenvalues of Φ are less than 1; it is unstable otherwise.

The step 4 is:

based on fixed mapping point theorem, the state variable vector x(t) satisfies x(t_(i))=x(t_(i)−T) for chatter-free cutting conditions; therefore, it has y_([0,T])=y_([−T,0]); state variables on all discrete points are calculated by:

$y^{*} = {y_{\lbrack{0,T}\rbrack} = {y_{\lbrack{{- T},0}\rbrack} = \left\{ \begin{matrix} {{\left( {P_{1} - P_{2} - Q} \right)^{- 1}z_{\lbrack{0,T}\rbrack}},} & {{{if}{❘{P_{1} - P_{2} - Q}❘}} \neq 0} \\ {{\left( {P_{1} - P_{2} - Q} \right)^{\dagger}z_{\lbrack{0,T}\rbrack}},} & {{{if}{❘{P_{1} - P_{2} - Q}❘}} = 0} \end{matrix} \right.}}$

where y* consists of modal displacement vector Γ(t_(i)) and modal velocity vector {dot over (Γ)}(t_(i)) at all discrete time nodes i=0,1, . . . ,2i_(m)+2 during one rotation period of cutting tool.

The step 5 comprises sub-steps of:

step 5.1: transform modal displacements into physical space, and calculate relative vibration displacements q(t_(i)) between cutting tool and workpiece:

q(t _(i))=q _(T)(t _(i))−q _(W)(t _(i))=P _(T)Γ_(T)(t _(i))−P _(W)Γ_(W)(t _(i))

where q_(T) (t_(i)) and q_(W)(t_(i)) are vibration displacements of cutting tool and workpiece, respectively;

step 5.2: extract vibration displacements in direction normal to milled surface of workpiece (defined as {right arrow over (Oy)} direction) from q(t_(i)), which constructs a new vector Q_(y)*:

Q _(y)*={[q _(y)(t ₁)]^(T),[q _(y)(t ₂)]^(T), . . . ,[q _(y)(t _(2i) _(m) ₊₂)]^(T)}

extract vibration displacements in feed direction (defined as {right arrow over (Ox)} direction), which constructs a new vector Q_(x)*:

Q _(x)*:={[q _(x)(t ₁)]^(T),[q _(x)(t ₂)]^(T), . . . ,[q _(x)(t _(2i) _(m) ₊₂)]^(T)}

step 5.3: based on kinematic synthesis of relative feed movement between cutting tool and workpiece, rotation of cutting tool and relative vibrations between cutting tool and workpiece, calculate motion trajectories in {right arrow over (Oy)} direction and {right arrow over (Ox)} direction for cutting element (j,k), respectively:

${\Theta_{y}^{\star}\left( {i,j,k} \right)} = {\underset{feed}{\underset{︸}{\frac{f \cdot t_{i}}{1000 \times 60}\sin\left( \varphi_{fy} \right)}} - \underset{rotation}{\underset{︸}{{R\left( {j,k} \right)}{\cos\left( {\varphi_{ref} - {\phi\left( {j,k} \right)} - \frac{z_{j}{\tan\left( \beta_{k} \right)}}{R} + {2\pi\frac{\Omega t_{i}}{60}}} \right)}}} + \underset{vibration}{\underset{︸}{Q_{y}^{*}\left( {i,j} \right)}}}$ ${\Theta_{x}^{\star}\left( {j,k,i} \right)} = {\underset{feed}{\underset{︸}{\frac{f \cdot t_{i}}{1000 \times 60}{\cos\left( \varphi_{fy} \right)}}} + \text{ }\underset{rotation}{\underset{︸}{{R\left( {j,k} \right)}{\sin\left( {\varphi_{ref} - {\phi\left( {j,k} \right)} - \frac{z_{j}{\tan\left( \beta_{k} \right)}}{R} + \frac{60t_{i}}{\Omega}} \right)}}} + \underset{vibrations}{\underset{︸}{Q_{x}^{*}\left( {j,i} \right)}}}$

where f is relative feed rate between cutting tool and workpiece in mm/min; φ_(fy) is angle between feed direction of cutting tool and normal direction of workpiece; R(j,k) is actual cutting radius of cutting element (j,k); φ_(ref) is reference angle used in the Generalized Runge-Kutta method; φ(j,k) is pitch angle between Tooth k and Tooth 1 for j-th discrete axial slice; β_(k) is helix angle with regard to Tooth k; z_(j) is axial height of j-th discrete axial slice; R is geometric radius of cutting tool; Ω is spindle rotation speed; Q_(y)*(j,i) and Q_(x)*(j,i) are vibration displacements of cutting element (j,k) in {right arrow over (Oy)} and {right arrow over (Ox)} directions, respectively.

The step 6 comprises sub-steps of:

step 6.1: select part of motion trajectories of cutting tool that are near milled surface of workpiece based on following equation, which construct motion trajectory points to be densified with interpolation, i.e., (Θ_(x_trim)*(j,k,i), ι_(y_trim)*(j,k,i)):

Θ_(y)*(j,k,i)≥max{Θ_(y)*(j,k,i)}−αR, down-milling

Θ_(y)*(j,k,i)≤min{Θ_(y)*(j,k,i)}+αR,up-milling

where α is a coefficient to adjust selecting range;

step 6.2: determine lower and upper limits of the selected trajectory points in {right arrow over (Ox)} direction, i.e., x_(min)=min {Θ_(x_trim)*(j,k,i)} and x_(max)=max {Θ_(x_trim)*(j,k,i)}; discretize interval [x_(min),x_(max)] equally with a step δx to obtain x coordinate set of interpolation points {x_(min), x_(min)+δx, x_(min)+2·δx, . . . , x_(max)}, where number of points is denoted by N_(s);

step 6.3: taking the selected points (Θ_(x_trim)*(j,k,i),Θ_(y_trim)*(j,k,i)) in step 6.1 as known points, use spline interpolation algorithm to determine y coordinates of interpolation points with regard to x coordinates {x_(min), x_(min)+δx,x_(min)+2·δx, . . . ,x_(max)}, which results in densified motion trajectories of cutting tool edges during one rotation period T of cutting tool, i.e., (x_(s)(l), y_(s)(l)),l=1, 2, . . . N_(s);

step 6.4: based on periodicity of chatter-free milling dynamic responses, i.e., x_(s)(l) and x_(s)(l)+n_(rev)T share same y coordinate y_(s)(l), obtain densified motion trajectory points of cutting tool during n_(rev), rotation periods of cutting tool, i.e., (x_(s_n)(l), y_(s_n)(l)), l=1,2, . . . N_(s_n), where N_(s_n) is number of interpolation points on n_(rev) rotation periods of cutting tool, and n_(rev)≥4;

step 6.5: select densified motion trajectory points of cutting tool edges on middle periods that satisfy x_(min)+T≤x_(s_n)(l)≤n_(rev)·x_(max)−T to construct a new densified point set (x_(s_n_trim)(l), y_(s_n_trim)(l)), l=1,2, . . . N_(s_n_trim), where N_(s_n_trim) is number of trimmed interpolation points;

step 6.6: sort the selected set (x_(s_n_trim)(l), y_(s_n_trim)(l)) in order of axial height to construct N_(a) sub-sets of densified points, i.e., (x_(s_n_trim),(j,l), y_(s_n_trim)(j,l)), j=1, 2, . . . ,N_(a);

step 6.7: at each axial height, compare values of y coordinates for all teeth that have same x coordinate; based on following equation, select densified points that are nearest to milled surface of workpiece to construct final surface topography points of workpiece (x_(surf)(j,l), y_(surf)(j,l)), l=1, 2, . . . N_(s_n_trim):

$\begin{matrix} {{{y_{surf}\left( {j,l} \right)} = {\max\limits_{k}\left\{ {y_{{s\_ n}{\_{trim}}}\left( {j,l} \right)} \right\}}},} & {{down} - {milling}} \end{matrix}$ $\begin{matrix} {{{y_{surf}\left( {j,l} \right)} = {\min\limits_{k}\left\{ {y_{{s\_ n}{\_{trim}}}\left( {j,l} \right)} \right\}}},} & {{up} - {milling}} \end{matrix}$

The step 7 is:

calculate surface location error (SLE) based on motion trajectories Θ_(y)*(i,j,k) in normal direction of workpiece; surface location error SLE(j,k) with regard to cutting element (j,k) is:

${SL{E\left( {j,k} \right)}} = \left\{ \begin{matrix} {{{\max\limits_{i}\left\{ {{\Theta_{y}^{*}\left( {j,k,i} \right)}❘{1 \leq i \leq {{2i_{m}} + 2}}} \right\}} - {D/2}},} & {{down} - {milling}} \\ {{{{- \min\limits_{i,k}}\left\{ {{\Theta_{y}^{*}\left( {j,k,i} \right)}❘{1 \leq i \leq {{2i_{m}} + 2}}} \right\}} - {D/2}},} & {{up}‐{milling}} \end{matrix} \right.$

surface location error with regard to axial height z_(j) is:

${SL{E(j)}} = \left\{ \begin{matrix} {{{\max\limits_{i,k}\left\{ {{{\Theta_{y}^{*}\left( {j,k,i} \right)}❘{1 \leq i \leq {{2i_{m}} + 2}}},{1 \leq k \leq N}} \right\}} - {D/2}},} & {{down} - {milling}} \\ {{{{- \min\limits_{i,k}}\left\{ {{{\Theta_{y}^{*}\left( {j,k,i} \right)}❘{1 \leq i \leq {{2i_{m}} + 2}}},{1 \leq k \leq N}} \right\}} - {D/2}},} & {{up}‐{milling}} \end{matrix} \right.$

where positive value mean overcut, and negative value mean undercut.

The step 8 is:

based on points (x_(surf)(j,l), y_(surf)(j,l)) that participate in generating surface topography of workpiece, surface roughness is calculated by:

${R{a(j)}} = {\frac{1}{N_{{s\_ n}{\_{trim}}}}{\sum\limits_{i^{\prime} = 1}^{N_{{s\_ n}{\_{trim}}}}{❘{{y_{surf}\left( {j,l} \right)} - {{\overset{¯}{y}}_{surf}\left( {j,l} \right)}}❘}}}$

By comparison with existing techniques, the invention has following advantages:

1. The invention discloses a method for simulating chatter-free milled surface topography. By comparison with initial value based time-domain simulation method, the invented method significantly improves the computational efficiency for simulating chatter-free milled surface topography;

2. The invention can simultaneously predict milling stability, surface location error and surface roughness, which provides theoretical and technical guidance for avoidance of milling chatter, improvement of cutting efficiency and assurance of machining quality.

BRIEF DESCRIPTION OF DRAWINGS

Other features, objects and advantages of the present invention will become apparent upon reading following detailed description of unrestricted examples, while referring to attached drawings in which:

FIGS. 1(a) to 1(d) are a flow diagram representing simulation procedure of milled surface topography of workpiece, where FIG. 1(a) represents motion trajectories of cutting tool edges, FIG. 1(b) represents densified motion trajectories of cutting tool using spline interpolation for each tooth during one rotation period of cutting tool, FIG. 1(c) represents densified motion trajectories of cutting tool during several rotation periods of cutting tool, and FIG. 1(d) represents the final surface topography of workpiece, which is determined by comparing motion trajectories of cutting tool edges for different teeth.

FIGS. 2(a) to 2(b) are comparison between experimental surface topography and simulated results; FIG. 2(a) is microphotograph of milled surface topography, and FIG. 2(b) is simulated milled surface topography.

FIG. 3 is comparison between simulated and measured surface location errors along tool axial height.

FIG. 4 is simulated surface roughness value along tool axial height.

EXAMPLES

The present invention will be described in detail below in conjunction with specific embodiments. The following examples will help those skilled in the art to further understand the present invention, but do not limit the present invention in any way. It should be pointed out that for those of ordinary skill in the art, several modifications and improvements can be made without departing from the concept of the present invention. These all belong to the protection scope of the present invention.

Refer to FIG. 1(a), FIG. 1(b), FIG. 1(c), FIG. 1(d), FIG. 2(a) and FIG. 2(b).

This example provides a method for simulating chatter-free milled surface topography, comprising steps of:

Step 1 comprises sub-steps of:

step 1.1: considering flexibility of both cutting tool and workpiece, variance of pitch angles and helix angles of cutting tool, and teeth runout, perform dynamics modeling on milling system and build a delay differential equation with multiple delays to represent the dynamic model in physical space:

$\begin{matrix} {{{M{\overset{¨}{q}(t)}} + {C{\overset{˙}{q}(t)}} + {K{q(t)}}} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{{K_{c}\left( {t,j,k} \right)}\left\lbrack {{q(t)} - {q\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack} + {F_{0}\left( {t,j,k} \right)}} \right\}}}} & (1) \end{matrix}$

where M, C and K are mass, damping and stiffness matrices, respectively; {umlaut over (q)}(t) {dot over (q)}(t) and q(t) are acceleration, velocity and displacement vectors at time t, respectively; K_(c)(t,j,k) is directional coefficient matrix for Tooth k at j-th axial slice at time t; F₀(t,j,k) is cutting force vector for Tooth k at j-th axial slice at time t, which is independent of chip regeneration; Σ_(p=k) _(v) ^(k)τ(j,p) is time delay with regard to cutting element (j,k), Σ is mathematical summation operator, p is a variable for summation, and k_(v) is tooth serial number to start delay calculation; N is number of teeth; and N_(a) is number of axial discretization of tool;

step 1.2: perform modal coordinate transformation on the dynamic model in step 1.1, which can be transformed from physical space into modal space; the modal coordinate transformation equation is:

q(t)=PΓ(t)  (2)

where P is the mode shape matrix, and Γ(t) is the modal displacement at time t;

dynamics model in modal space which is represented with delay differential equation with multiple delays turns out to be:

$\begin{matrix} {{{M_{\Gamma}{\overset{¨}{\Gamma}(t)}} + {C_{\Gamma}{\overset{.}{\Gamma}(t)}} + {K_{\Gamma}{\Gamma(t)}}} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{P^{T}{K_{c}\left( {t,j,k} \right)}{P\left\lbrack {{\Gamma(t)} - {\Gamma\left( {t - {\sum_{p = k_{v}}^{k}\ {\tau\left( {j,p} \right)}}} \right)}} \right\rbrack}} + {P^{T}{F_{0}\left( {t,j,k} \right)}}} \right\}}}} & (3) \end{matrix}$

where M_(Γ), C_(Γ) and K_(Γ) are modal mass, modal damping and modal stiffness matrices, respectively; {umlaut over (Γ)}(t), {dot over (Γ)}(t) and Γ(t) are modal acceleration, modal velocity and modal displacement vectors, respectively;

for milling system with flexible cutting tool and flexible workpiece, the dynamic model becomes:

$\begin{matrix} {{{\begin{bmatrix} M_{\Gamma;T} & \\  & M_{\Gamma;W} \end{bmatrix}\begin{Bmatrix} {{\overset{¨}{\Gamma}}_{T}(t)} \\ {{\overset{¨}{\Gamma}}_{W}(t)} \end{Bmatrix}} + {\begin{bmatrix} C_{\Gamma;T} & \\  & C_{\Gamma;W} \end{bmatrix}\begin{Bmatrix} {{\overset{.}{\Gamma}}_{T}(t)} \\ {{\overset{.}{\Gamma}}_{W}(t)} \end{Bmatrix}} + \begin{bmatrix} K_{\Gamma;T} & \\  & K_{\Gamma;W} \end{bmatrix}}\text{ }{\begin{Bmatrix} {\Gamma_{T}(t)} \\ {\Gamma_{W}(t)} \end{Bmatrix} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{\begin{bmatrix} P_{T}^{T} \\ {- P_{W}^{T}} \end{bmatrix}{{K_{c}\left( {t,j,k} \right)}\begin{bmatrix} P_{T} & {- P_{W}} \end{bmatrix}}\text{ }\begin{Bmatrix} {{\Gamma_{T}(t)} - {\Gamma_{T}\left( {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}} \right)}} \\ {{\Gamma_{W}(t)} - {\Gamma_{W}\left( {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}} \right)}} \end{Bmatrix}} + {\begin{bmatrix} P_{T}^{T} \\ {- P_{W}^{T}} \end{bmatrix}{F_{0}\left( {t,j,k} \right)}}} \right\}}}}} & (4) \end{matrix}$

where M_(Γ,T), C_(Γ,T) and K_(Γ,T) are modal mass, modal damping and modal stiffness matrices of cutting tool; M_(Γ;W), C_(Γ;W) and K_(Γ;W) are modal mass, modal damping and modal stiffness matrices of workpiece; {umlaut over (Γ)}_(T)(t), {dot over (Γ)}_(T)(t) and Γ_(T)(t) are modal acceleration, modal velocity and modal displacement vectors of cutting tool; {umlaut over (Γ)}_(W)(t), {dot over (Γ)}_(W)(t) and Γ_(W)(t) are modal acceleration, modal velocity and modal displacement vectors of workpiece; P_(T) and P_(W) are mode shape matrices of cutting tool and workpiece, respectively; P_(T) ^(T) and P_(W) ^(T) are transposed matrices of P_(T) and P_(W), respectively;

for milling system with flexible cutting tool and rigid workpiece, the dynamic model becomes:

$\begin{matrix} {{{M_{\Gamma;T}{{\overset{¨}{\Gamma}}_{T}(t)}} + {C_{\Gamma;T}{{\overset{˙}{\Gamma}}_{T}(t)}} + {K_{\Gamma;T}{\Gamma_{T}(t)}}} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{P_{T}^{T}{K_{c}\left( {t,j,k} \right)}\left\{ {P_{T}\left\lbrack {{\Gamma_{T}(t)} - {\Gamma_{T}\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack} \right\}} + {P_{T}^{T}{F_{0}\left( {t,j,k} \right)}}} \right\}}}} & (5) \end{matrix}$

for milling system with rigid cutting tool and flexible workpiece, the dynamic model becomes:

$\begin{matrix} {{{M_{\Gamma;W}{{\overset{¨}{\Gamma}}_{W}(t)}} + {C_{\Gamma;W}{{\overset{˙}{\Gamma}}_{W}(t)}} + {K_{\Gamma;W}{\Gamma_{W}(t)}}} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{P_{W}^{T}{K_{c}\left( {t,j,k} \right)}\left\{ {P_{W}\left\lbrack {{\Gamma_{W}(t)} - {\Gamma_{W}\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack} \right\}} - {P_{W}^{T}{F_{0}\left( {t,j,k} \right)}}} \right\}}}} & (6) \end{matrix}$

Step 2 comprises sub-steps of:

step 2.1: perform state space transformation on the dynamic model in modal space by introducing state variable vector x(t)=[(Γ(t))^(T), ({dot over (Γ)}(t))^(T)]^(T), and delay differential equation with multiple delays in state space becomes:

$\begin{matrix} {{\overset{˙}{x}(t)} = {{A{x(t)}} + {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{{B\left( {t,j,k} \right)}\left\lbrack {{x(t)} - {x\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack} + {D\left( {t,j,k} \right)}} \right\}}}}} & (7) \end{matrix}$

where x(t) denotes state variable vector, {dot over (x)}(t) denotes first derivative of x(t) with regard to time t,

${A = \begin{bmatrix} 0 & I \\ {{- M_{\Gamma}^{- 1}}K_{\Gamma}} & {{- M_{\Gamma}^{- 1}}C_{\Gamma}} \end{bmatrix}}\ ,$ ${{B\left( {t,j,k} \right)} = \begin{bmatrix} 0 & 0 \\ {M_{\Gamma}^{- 1}P^{T}{K_{c}\left( {t,j,k} \right)}P} & 0 \end{bmatrix}},$ ${{D\left( {t,j,k} \right)} = \begin{bmatrix} 0 \\ {M_{\Gamma}^{- 1}P^{T}{F_{0}\left( {t,j,k} \right)}} \end{bmatrix}},$

and I is unit matrix;

step 2.2: divide one rotation period of cutting tool T into 2i_(m)+2 equal intervals with a series of discrete nodes {t₀,t₁, . . . ,t_(2i) _(m) ₊₂}; analytical expression of dynamic responses at interval [t_(i),t] with regard to the delay differential equation with multiple delays in step 2.1 is:

$\begin{matrix} {{x(t)} = {{e^{A({t - t_{i}})}{x\left( t_{i} \right)}} + \text{ }{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}{\int_{t_{i}}^{t}{\left\{ {{e^{A({t - \xi})}{{B\left( {\xi,j,k} \right)}\left\lbrack {{x(\xi)} - {x\left( {\xi - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack}} + \text{ }{e^{A({t - \xi})}{D\left( {\xi,j,k} \right)}}} \right\} d{\xi.}}}}}}} & (8) \end{matrix}$

for brevity, define x_(i):=x(t_(i)), E(t,t_(i)):=e^(A(t−t) ^(i) ⁾x(t_(i)), H_(j,k)(t,ξ,x(ξ)):=e^(A(t−ξ))B(ξ, j,k)[x(ξ)−x(ξ−τ(j,k, p))] and J_(j,k)(t,ξ):=e^(A(t−ξ))D(ξ,j, j,k);

on sub-interval [t_(2i),t_(2i+2)], x(t) is approximated using Simpson's rule:

$\begin{matrix} {x_{{2i} + 2} = {{E\left( {t_{{2i} + 2},t_{2i}} \right)} + {\frac{h}{3}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{H_{j,k}\left( {t_{{2i} + 2},t_{2i},x_{2i}} \right)} + \text{ }{4{H_{j,k}\left( {t_{{2i} + 2},t_{{2i} + 1},x_{{2i} + 1}} \right)}} + {H_{j,k}\left( {t_{{2i} + 2},t_{{2i} + 2},x_{{2i} + 2}} \right)} + {J_{j,k}\left( {t_{{2i} + 2},t_{2i}} \right)} + {4{J_{j,k}\left( {t_{{2i} + 2},t_{{2i} + 1}} \right)}} + {J_{j,k}\left( {t_{{2i} + 2},t_{{2i} + 2}} \right)}} \right\}}}}}} & (9) \end{matrix}$

on sub-interval [t_(2i),t_(2i+1)], x(t) is approximated using fourth-order Runge-Kutta algorithm:

$\begin{matrix} {x_{{2i} + 1} = {{E\left( {t_{{2i} + 1},t_{2i}} \right)} + {\frac{h}{6}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{H_{j,k}\left( {t_{{2i} + 1},t_{2i},x_{2i}} \right)} + {2{H_{j,k}\left( {t_{{2i} + 1},t_{{2i} + {1/2}},x_{{2i} + {1/2}}} \right)}} + {2{H_{j,k}\left( {t_{{2i} + 1},t_{{2i} + {1/2}},x_{{2i} + {1/2}}} \right)}} + {H_{j,k}\left( {t_{{2i} + 1},t_{{2i} + 1},x_{{2i} + 1}} \right)} + {J_{j,k}\left( {t_{{2i} + 1},t_{2i}} \right)} + {2{J_{j,k}\left( {t_{{2i} + 1},t_{{2i} + {1/2}}} \right)}} + {2{J_{j,k}\left( {t_{{2i} + 1},t_{{2i} + {1/2}}} \right)}} + {J_{j,k}\left( {t_{{2i} + 1},t_{{2i} + 1}} \right)}} \right\}}}}}} & (10) \end{matrix}$

where state variable vector at middle points t_(2i+1/2) is approximated with Lagrange formula:

x _(2i+1/2)=⅜x _(2i)−¾x _(2i−1)+⅛x _(2i+2)  (11)

the analytical expression of dynamic responses on sub-interval [t_(2i),t_(2i+1)] is re-organized into:

$\begin{matrix} \begin{matrix} {{{F_{2i}x_{2i}} + {F_{{2i} + 1}x_{{2i} + 1}} + {F_{{2i} + 2}x_{{2i} + 2}}} =} \\ {{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left( {{F_{{2i} - {\tau({j,k})}}x_{{2i} - {\tau({j,k})}}} + {F_{{2i} + 1 - {\tau({j,k})}}x_{{2i} + 1 - {\tau({j,k})}}} + {F_{{2i} + 2 - {\tau({j,k})}}x_{{2i} + 2 - {\tau({j,k})}}}} \right)}} + {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left( {S_{2i} + S_{{2i} + {1/2}} + S_{{2i} + 1}} \right)}}} \\ {where} \\ {{F_{2i} = {{- e^{Ah}} - {\frac{h}{6}e^{Ah}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{2i},j,k}}}} - {{\frac{3}{8} \cdot \frac{4h}{6}}e^{Ah/2}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + {1/2}},j,k}}}}}},} \\ {{F_{{2i} + 1} = {I - {{\frac{3}{4} \cdot \frac{4h}{6}}e^{Ah/2}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + {1/2}},j,k}}}} - {\frac{h}{6}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + 1},j,k}}}}}},{F_{{2i} + 2} = {{\frac{1}{8} \cdot \frac{4h}{6}}e^{Ah/2}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + {1/2}},j,k}}}}},} \\ {{F_{{2i} - {\tau({j,k,p})}} = {{- \frac{h}{6}e^{Ah}B_{{2i},j,k}} - {{\frac{3}{8} \cdot \frac{4h}{6}}e^{Ah/2}B_{{{2i} + {1/2}},j,k}}}},{F_{{2i} + 1 - {\tau({j,k,p})}} = {{- {\frac{3}{4} \cdot \frac{4h}{6}}e^{Ah/2}B_{{{2i} + {1/2}},j,k}} - {\frac{h}{6}B_{{{2i} + 1},j,k}}}},} \\ {{F_{{2i} + 2 - {\tau({j,k,p})}} = {{\frac{1}{8} \cdot \frac{4h}{6}}e^{Ah/2}B_{{{2i} + {1/2}},j,k}}},{S_{2i} = {\frac{h}{6}e^{A \cdot h}D_{{2i},j,k}}},{S_{{2i} + {1/2}} = {{\frac{h}{6} \cdot 4}e^{{A \cdot h}/2}D_{{{2i} + {1/2}},j,k}}},{S_{{2i} + 1} = {\frac{h}{6}D_{{{2i} + 1},j,k}}},} \end{matrix} & (12) \end{matrix}$

and h is discretization step;

the analytical expression of dynamic responses on sub-interval [t_(2i),t_(2i+2)] is re-organized into:

$\begin{matrix} \begin{matrix} {{{G_{2i}x_{2i}} + {G_{{2i} + 1}x_{{2i} + 1}} + {G_{{2i} + 2}x_{{2i} + 2}}} =} \\ {{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left( {{G_{{2i} - {\tau({j,k})}}x_{{2i} - {\tau({j,k})}}} + {G_{{2i} + 1 - {\tau({j,k})}}x_{{2i} + 1 - {\tau({j,k})}}} + {G_{{2i} + 2 - {\tau({j,k})}}x_{{2i} + 2 - {\tau({j,k})}}}} \right)}} + {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left( {T_{2i} + T_{{2i} + 1} + T_{{2i} + 2}} \right)}}} \\ {where} \\ {{G_{2i} = {{- e^{2Ah}} - {\frac{h}{3}e^{2Ah}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{2i},j,k}}}}}},{G_{{2i} + 1} = {- \frac{4h}{3}e^{Ah}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + 1},j,k}}}}},} \\ {{G_{{2i} + 2} = {I - {\frac{h}{3}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + 2},j,k}}}}}},{G_{{2i} - {\tau({j,k,p})}} = {- \frac{h}{3}e^{2Ah}B_{{2i},j,k}}},{G_{{2i} + 1 - {\tau({j,k,p})}} = {- \frac{4h}{3}e^{Ah}B_{{{2i} + 1},j,k}}},} \\ {{G_{{2i} + 2 - {\tau({j,k,p})}} = {- \frac{h}{3}B_{{{2i} + 2},j,k}}},{T_{2i} = {\frac{h}{3}e^{{A \cdot 2}h}D_{{2i},j,k}}},{{T_{{2i} + 1} = {{{\frac{h}{3} \cdot 4}e^{A \cdot h}D_{{{2i} + 1},j,k}{and}T_{{2i} + 2}} = {\frac{h}{3}D_{{{2i} + 2},j,k}}}};}} \end{matrix} & (13) \end{matrix}$

step 2.3: based on the formulae in step 2.2, discrete mapping relationship between two adjacent rotation periods of cutting tool [−T,0] and [0,T] is established:

P ₁ y _([0,T]) =Qy _([−T,0]) +P ₂ y _([0,T]) +z _([0,T])  (14)

where

$\begin{matrix} {{y_{\lbrack{0,T}\rbrack} = \begin{bmatrix} x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\  \vdots \\ x_{{2i_{m}} + 1} \\ x_{{2i_{m}} + 2} \end{bmatrix}},{y_{{\lbrack{{- T},0}\rbrack} =}\begin{bmatrix} x_{0 - T} \\ x_{1 - T} \\ x_{2 - T} \\ x_{3 - T} \\ x_{4 - T} \\  \vdots \\ x_{{2i_{m}} + 1 - T} \\ x_{{2i_{m}} + 2 - T} \end{bmatrix}},{z_{\lbrack{0,T}\rbrack} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\begin{bmatrix} 0 \\ {S_{0} + S_{1/2} + S_{1}} \\ {T_{0} + T_{1} + T_{2}} \\ {S_{2} + S_{2 + {1/2}} + S_{3}} \\ {T_{2} + T_{3} + T_{4}} \\  \vdots \\ {S_{2i_{m}} + S_{{2i_{m}} + {1/2}} + S_{{2i_{m}} + 1}} \\ {T_{2i_{m}} + T_{{2i_{m}} + 1} + T_{{2i_{m}} + 2}} \end{bmatrix}}}},} \\ {{P_{1} = \begin{bmatrix} I & 0 & 0 & 0 & 0 & \ldots & 0 & 0 & 0 \\ F_{0} & F_{1} & F_{2} & 0 & 0 & \ldots & 0 & 0 & 0 \\ G_{0} & G_{1} & G_{2} & 0 & 0 & \ldots & 0 & 0 & 0 \\ 0 & 0 & F_{2} & F_{3} & F_{4} & \ldots & 0 & 0 & 0 \\ 0 & 0 & G_{2} & G_{3} & G_{4} & \ldots & 0 & 0 & 0 \\  \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & 0 & \ldots & F_{2i_{m}} & F_{{2i_{m}} + 1} & F_{{2i_{m}} + 2} \\ 0 & 0 & 0 & 0 & 0 & \ldots & G_{2i_{m}} & G_{{2i_{m}} + 1} & G_{{2i_{m}} + 2} \end{bmatrix}},} \\ {Q = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\begin{bmatrix} 0 & \ldots & 0 & 0 & 0 & \ldots & 0 & 0 & {I/{NN}_{a}} \\ 0 & \ldots & F_{- {\tau({j,k,p})}} & F_{1 - {\tau({j,k,p})}} & F_{2 - {\tau({j,k,p})}} & \ldots & 0 & 0 & 0 \\ 0 & \ldots & G_{- {\tau({j,k,p})}} & G_{1 - {\tau({j,k,p})}} & G_{2 - {\tau({j,k,p})}} & \ddots & \vdots & 0 & 0 \\  \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & \ldots & 0 & 0 & 0 & \ldots & F_{{2i} - {\tau({j,k,p})}} & F_{{2i} + 1 - {\tau({j,k,p})}} & F_{{2i} + 2 - {\tau({j,k,p})}} \\ 0 & \ldots & 0 & 0 & 0 & \ldots & G_{{2i} - {\tau({j,k,p})}} & G_{{2i} + 1 - {\tau({j,k,p})}} & G_{{2i} + 2 - {\tau({j,k,p})}} \\  \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & \ldots & 0 & 0 & 0 & \ldots & 0 & 0 & 0 \end{bmatrix}}}} \\ {and} \\ {P_{2} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}{\begin{bmatrix} 0 & 0 & 0 & \ldots & 0 & 0 & 0 & \ldots & 0 \\  \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ldots & 0 \\ F_{{2i} + 2 - {\tau({j,k,p})}} & F_{{2i} + 3 - {\tau({j,k,p})}} & F_{{2i} + 4 - {\tau({j,k,p})}} & \ldots & 0 & 0 & 0 & \ldots & 0 \\ G_{{2i} + 2 - {\tau({j,k,p})}} & G_{{2i} + 3 - {\tau({j,k,p})}} & G_{{2i} + 4 - {\tau({j,k,p})}} & \ddots & 0 & 0 & 0 & \ldots & 0 \\ 0 & 0 & \vdots & \ddots & 0 & 0 & 0 & \ldots & 0 \\  \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \ldots & F_{{2i_{m}} - {\tau({j,k,p})}} & F_{{2i_{m}} + 1 - {\tau({j,k,p})}} & F_{{2i_{m}} + 2 - {\tau({j,k,p})}} & \ldots & 0 \\ 0 & 0 & 0 & \ldots & G_{{2i_{m}} - {\tau({j,k,p})}} & G_{{2i_{m}} + 1 - {\tau({j,k,p})}} & G_{{2i_{m}} + 2 - {\tau({j,k,p})}} & \ldots & 0 \end{bmatrix}.}}}} \end{matrix}$

Step 3 is:

construct transition matrix Φ of the milling system:

$\begin{matrix} {\Phi = \left\{ \begin{matrix} {{\left( {P_{1} - P_{2}} \right)^{- 1}Q},{{{if}{❘{P_{1} - P_{2}}❘}} \neq 0}} \\ {{\left( {P_{1} - P_{2}} \right)^{\dagger}Q},{{{if}{❘{P_{1} - P_{2}}❘}} = 0}} \end{matrix} \right.} & (15) \end{matrix}$

where |P₁−P₂| denotes determinant of matrix (P₁−P₂); † denotes Moore-Penrose pseudoinverse of matrix (P₁−P₂); based on Floquet theory, stability of milling system depends on eigenvalue of Φ: it is stable if all modules of eigenvalues of Φ are less than 1; it is unstable otherwise.

Step 4 is:

based on fixed mapping point theorem, the state variable vector x(t) satisfies x(t_(i))=x(t_(i)−T) for chatter-free cutting conditions; therefore, it has y_([0,T])=y_([−T,0]); state variables on all discrete points are calculated by:

$\begin{matrix} {y^{*} = {y_{\lbrack{0,T}\rbrack} = {y_{\lbrack{{- T},0}\rbrack} = \left\{ \begin{matrix} {{\left( {P_{1} - P_{2} - Q} \right)^{- 1}z_{\lbrack{0,T}\rbrack}},{{{if}{❘{P_{1} - P_{2} - Q}❘}} \neq 0}} \\ {{\left( {P_{1} - P_{2} - Q} \right)^{\dagger}z_{\lbrack{0,T}\rbrack}},{{{if}{❘{P_{1} - P_{2} - Q}❘}} = 0}} \end{matrix} \right.}}} & (16) \end{matrix}$

where y* consists of modal displacement vector Γ(t_(i)) and modal velocity vector {dot over (Γ)}(t_(i)) at all discrete time nodes i=0,1, . . . ,2i_(m)30 2 during one rotation period of cutting tool.

Step 5 comprises sub-steps of:

step 5.1: transform modal displacements into physical space, and calculate relative vibration displacements q(t_(i)) between cutting tool and workpiece:

q(t _(i))=q _(T)(t _(i))−q _(W)(t _(i))=P _(T)Γ_(T)(t _(i))−P _(W)Γ_(W)(t _(i))  (17)

where q_(T)(t_(i)) and q_(W)(t_(i)) are vibration displacements of cutting tool and workpiece, respectively;

step 5.2: extract vibration displacements in direction normal to milled surface of workpiece (defined as {right arrow over (Oy)} direction) from q(t_(i)), which constructs a new vector Q_(y)*:

Q _(y)*={[q _(y)(t ₁)]^(T),[q _(y)(t ₂)]^(T), . . . ,[q _(y)(t _(2i) _(m) ₊₂)]^(T)}  (18)

extract vibration displacements in feed direction (defined as {right arrow over (Ox)} direction), which constructs a new vector Q_(x)*:

Q _(x)*={[q _(x)(t ₁)]^(T),[q _(x)(t ₂)]^(T), . . . ,[q _(x)(t _(2i) _(m) ₊₂)]^(T)}  (19)

step 5.3: based on kinematic synthesis of relative feed movement between cutting tool and workpiece, rotation of cutting tool and relative vibrations between cutting tool and workpiece, calculate motion trajectories in {right arrow over (Oy)} direction and {right arrow over (Ox)} direction for cutting element (j,k), respectively, as shown in FIG. 1 (a):

$\begin{matrix} {{\Theta_{y}^{*}\left( {i,j,k} \right)} = {\underset{feed}{\underset{︸}{\frac{f \cdot t_{i}}{1000 \times 60}\sin\left( \varphi_{fy} \right)}} - \underset{rotation}{\underset{︸}{\begin{matrix} {{R\left( {j,k} \right)}{\cos\left( {\varphi_{ref} - {\phi\left( {j,k} \right)} - \frac{z_{j}{\tan\left( \beta_{k} \right)}}{R}} \right.}} & \left. {{+ 2}\pi\frac{\Omega t_{i}}{60}} \right) \end{matrix}}} + \underset{vibration}{\underset{︸}{\begin{matrix} Q_{y}^{*} & \left( {i,j} \right) \end{matrix}}}}} & (20) \end{matrix}$ $\begin{matrix} {{\Theta_{x}^{*}\left( {j,k,i} \right)} = {\underset{feed}{\underset{︸}{\begin{matrix} \frac{f \cdot t_{i}}{1000 \times 60} & {\cos\left( \varphi_{fy} \right)} \end{matrix}}} + \underset{rotation}{\underset{︸}{{R\left( {j,k} \right)}{\sin\left( {\varphi_{ref} + {\phi\left( {j,k} \right)} - \frac{z_{j}{\tan\left( \beta_{k} \right)}}{R} + \frac{60t_{i}}{\Omega}} \right)}}} + \underset{vibrations}{\underset{︸}{\begin{matrix} Q_{x}^{*} & \left( {j,i} \right) \end{matrix}}}}} & (21) \end{matrix}$

where f is relative feed rate between cutting tool and workpiece in mm/min; φ_(fy) is angle between feed direction of cutting tool and normal direction of workpiece; R(j,k) is actual cutting radius of cutting element (j,k); q); φ_(ref) is reference angle used in the Generalized Runge-Kutta method; ϕ(j,k) is pitch angle between Tooth k and Tooth 1 for j-th discrete axial slice; β_(k) is helix angle with regard to Tooth k; z_(j) is axial height of j-th discrete axial slice; R is geometric radius of cutting tool; Ω is spindle rotation speed; Q_(y)*(j,i) and Q_(x)*(j,i) are vibration displacements of cutting element (j,k) in {right arrow over (Oy)} and {right arrow over (Ox)} directions, respectively.

Step 6 comprises sub-steps of:

step 6.1: select part of motion trajectories of cutting tool that are near milled surface of workpiece based on following equation, which construct motion trajectory points to be densified with interpolation, i.e., (Θ_(x_trim)*(j,k,i), Θ_(y_trim)*(j,k,i)):

Θ_(y)*(j,k,i)≥max{Θ_(y)*(j,k,i)}−αR, down-milling

Θ_(y)*(j,k,i)≤min{Θ_(y)*(j,k,i)}+αR, up-milling

where α is a coefficient to adjust selecting range;

step 6.2: determine lower and upper limits of the selected trajectory points in {right arrow over (Ox)} direction, i.e., x_(min)=min {Θ_(x_trim)*(j,k,i)} and x_(max)=max {Θ_(x_trim)*(j,k,i)}; discretize interval [x_(min),x_(max)] equally with a step δx to obtain x coordinate set of interpolation points {x_(min),x_(min)+δx, x_(min)+2·δx, . . . ,x_(max)}, where number of points is denoted by N_(s);

step 6.3: taking the selected points (Θ_(x_trim)*(j,k,i), Θ_(y_trim)*(j,k,i)) in step 6.1 as known points, use spline interpolation algorithm to determine y coordinates of interpolation points with regard to x coordinates {x_(min),x_(min)+δx,x_(min)+2·δx, . . . ,x_(max)}, which results in densified motion trajectories of cutting tool edges during one rotation period T of cutting tool, i.e., (x_(s)(l), y_(s) (l)),l=1, 2, . . . N_(s), as shown in FIG. 1 (b);

step 6.4: based on periodicity of chatter-free milling dynamic responses, i.e., x_(s)(l) and x_(s)(l)+n_(rev)T share same y coordinate y_(s)(l), obtain densified motion trajectory points of cutting tool during n_(rev) rotation periods of cutting tool, i.e., (x_(s_n_trim)(l), y_(s_n_trim)(l)), l=1, 2, . . . N_(s_n_trim), where N_(s_n_trim) is number of interpolation points on n_(rev) rotation periods of cutting tool, and n_(rev)≥4, as shown in FIG. 1 (c);

step 6.5: select densified motion trajectory points of cutting tool edges on middle periods that satisfy x_(min), +T≤x_(s_n)(l)≤n_(rev)·x_(max)−T to construct a new densified point set (x_(s_n_trim)(l), y_(s_n_trim)(l)), l=1, 2, N_(s_n_trim), where N_(s_n_trim) is number of trimmed interpolation points;

step 6.6: sort the selected set (x_(s_n_trim) (l), y_(s_n_trim)(l)) in order of axial height to construct N_(a) sub-sets of densified points, i.e., (x_(s_n_trim)(j,l), y_(s_n_trim)(j,l), j=1, 2, . . . ,N_(a);

step 6.7: at each axial height, compare values of y coordinates for all teeth that have same x coordinate; based on following equation, select densified points that are nearest to milled surface of workpiece to construct final surface topography points of workpiece (x_(surf)(j,l), y_(surf)(j,l)), l=1,2, . . . N_(s_n_trim), shown in FIG. 1 (d) and FIG. 2 (b):

$\begin{matrix} {{{y_{surf}\left( {j,l} \right)} = {\max\limits_{k}\left\{ {y_{{s\_ n}{\_{trim}}}\left( {j,l} \right)} \right\}}},{{down} - {milling}}} \\ {{{y_{surf}\left( {j,l} \right)} = {\min\limits_{k}\left\{ {y_{{s\_ n}{\_{trim}}}\left( {j,l} \right)} \right\}}},{{up} - {milling}}} \end{matrix}$

Step 7 is:

calculate surface location error (SLE) based on motion trajectories Θ_(y)*(i,j,k) in normal direction of workpiece; surface location error SLE(j,k) with regard to cutting element (j,k) is:

$\begin{matrix} {{{SLE}\left( {j,k} \right)} = \left\{ \begin{matrix} {{{\max\limits_{i}\left\{ {{\Theta_{y}^{*}\left( {j,k,i} \right)}❘{1 \leq i \leq {{2i_{m}} + 2}}} \right\}} - {D/2}},{{down} - {milling}}} \\ {{{- \min\limits_{i}\left\{ {{\Theta_{y}^{*}\left( {j,k,i} \right)}❘{1 \leq i \leq {{2i_{m}} + 2}}} \right\}} - {D/2}},{{up} - {milling}}} \end{matrix} \right.} & (22) \end{matrix}$

surface location error with regard to axial height z_(j) is:

$\begin{matrix} {{{SLE}(j)} = \left\{ \begin{matrix} {{{\max\limits_{i,k}\left\{ {{{\Theta_{y}^{*}\left( {j,k,i} \right)}❘{1 \leq i \leq {{2i_{m}} + 2}}},{1 \leq k \leq N}} \right\}} - {D/2}},{{down} - {milling}}} \\ {{{- \min\limits_{i,k}\left\{ {{{\Theta_{y}^{*}\left( {j,k,i} \right)}❘{1 \leq i \leq {{2i_{m}} + 2}}},{1 \leq k \leq N}} \right\}} - {D/2}},{{up} - {milling}}} \end{matrix} \right.} & (23) \end{matrix}$

where positive value mean overcut, and negative value mean undercut, as shown in FIG. 3.

Step 8 is:

based on points (x_(surf) (j,l), y_(surf)(j,l)) that participate in generating surface topography of workpiece, surface roughness is calculated by:

$\begin{matrix} {{{Ra}(j)} = {\frac{1}{N_{{s\_ n}{\_ trim}}}{\sum\limits_{i^{\prime} = 1}^{N_{{s\_ n}{\_{trim}}}}{❘{{y_{surf}\left( {j,l} \right)} - {{\overset{\_}{y}}_{surf}\left( {j,l} \right)}}❘}}}} & (24) \end{matrix}$

An experimental example is presented to further illustrate the invented method. The geometric parameters of cutting tool are: diameter D=12 mm, number of teeth N=4, pitch angles 80°-100°-80°-100°, and helix angles 45°-45°-45°-45°. The cutting parameters are: spindle rotation speed Ω=2500 rpm, radial depth of cut a_(r)=0.5 mm, axial depth of cut a_(p)=5 mm, and feed per revolution f_(rev)=0.2 mm/rev. The system dynamic parameters are: natural frequency f_(n)=189.1 Hz, damping ratio

=0.0047, and stiffness k=3.01×10⁶ N/m. The runout parameters are: offset value ρ=2.5 μm, and orientation angle λ=0.1°.

With the provided parameters, milled surface topography can be simulated by following the procedures in steps 1˜7, as illustrated in FIG. 2. According to microphotography in FIG. 2(b), the simulated surface topography matches milled surface topography well. Simulated and measured surface location errors are illustrated in FIG. 3, where positive value means overcut and negative value means undercut. According to FIG. 3, the simulated surface location errors locate in the range between −43.5 μm and 32.6 μm. Selecting six points along axial height of cutting tool, measured surface location errors using Coordinate Measurement Machine are −41.6 μm, −29.0 μm, −18.7 μm, −11.1 μm, −6.2 μm and 0.3 μm, which validate the simulated results. As shown in FIG. 4, the predicted surface roughness locates in the range between 0.0679 μm and 0.2449 μm. Selecting four points along axial height of cutting tool, measured surface roughness values are R_(a)=0.2190 μm, R_(a)=0.2865 μm, R_(a)=0.2395 μm and R_(a)=0.2665 μm, which validate the simulated results.

This invention can be embodied in other forms without departing from the spirit or essential attributes thereof and, accordingly, reference should be had to the claims rather than the foregoing specification as indicating the scope of the invention. 

1. A method for simulating chatter-free milled surface topography, comprising steps of: step 1: performing dynamics modeling on milling system and building a delay differential equation with multiple delays to represent dynamic model; step 2: constructing state transition matrix between two adjacent cutting tool rotation periods by extending Generalized Runge-Kutta method; step 3: predicting stable milling parameter zones based on Floquet theory; step 4: calculating vibration displacements on discretized time nodes during one cutting tool rotation period by means of fixed mapping point theorem; step 5: constructing motion trajectories of cutting tool edges in normal and feed directions with regard to milled surface of workpiece; step 6: predicting surface topography of workpiece by using spline interpolation to densify motion trajectories of cutting tool edges that participate in generating final surface of workpiece; step 7: calculating surface location error of milled surface based on motion trajectories of cutting tool edges in normal direction with regard to milled surface of workpiece; step 8: calculating surface roughness based on the simulated surface topography of workpiece.
 2. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 1 comprises sub-steps of: step 1.1: considering flexibility of both cutting tool and workpiece, variance of pitch angles and helix angles of cutting tool, and teeth runout, performing dynamics modeling on milling system and building a delay differential equation with multiple delays to represent the dynamic model in physical space: ${{M{\overset{¨}{q}(t)}} + {C{\overset{.}{q}(t)}} + {{Kq}(t)}} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{{K_{c}\left( {t,j,k} \right)}\left\lbrack {{q(t)} - {q\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack} + {F_{0}\left( {t,j,k} \right)}} \right\}}}$ where M C and K are mass, damping and stiffness matrices, respectively; {umlaut over (q)}(t), {dot over (q)}(t) and q(t) are acceleration, velocity and displacement vectors at time t, respectively; K_(c)(t,j,k) is directional coefficient matrix for Tooth k at j-th axial slice at time t; F₀(t,j,k) is cutting force vector for Tooth k at j-th axial slice at time t, which is independent of chip regeneration; Σ_(p=k) ^(k)τ(j,p) is time delay with regard to cutting element (j,k) is mathematical summation operator, p is a variable for summation, and is tooth serial number to start delay calculation; N is number of teeth; and N_(a) is number of axial discretization of tool; step 1.2: performing modal coordinate transformation on the dynamic model in step 1.1, which can be transformed from physical space into modal space; the modal coordinate transformation equation is: q(t)=PΓ(t) where P is mode shape matrix, and Γ(t) is modal displacement vector at time t; dynamics model in modal space which is represented with delay differential equation with multiple delays turns out to be: ${{M_{\Gamma}{\overset{¨}{\Gamma}(t)}} + {C_{\Gamma}{\overset{.}{\Gamma}(t)}} + {K_{\Gamma}{\Gamma(t)}}} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{P^{T}{K_{c}\left( {t,j,k} \right)}{P\left\lbrack {{\Gamma(t)} - {\Gamma\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack}} + {P^{T}{F_{0}\left( {t,j,k} \right)}}} \right\}}}$ where M_(Γ), C_(Γ) and K_(Γ) are modal mass, modal damping and modal stiffness matrices, respectively; {umlaut over (Γ)}(t), {dot over (Γ)}(t) and Γ(t) are modal acceleration, modal velocity and modal displacement vectors, respectively; for milling system with flexible cutting tool and flexible workpiece, the dynamic model becomes: $\begin{matrix} {{{\begin{bmatrix} M_{\Gamma;T} & \text{ } \\ \text{ } & M_{\Gamma;W} \end{bmatrix}\begin{Bmatrix} {{\overset{¨}{\Gamma}}_{T}(t)} \\ {{\overset{¨}{\Gamma}}_{W}(t)} \end{Bmatrix}} + {\begin{bmatrix} C_{\Gamma;T} & \text{ } \\ \text{ } & C_{\Gamma;W} \end{bmatrix}\begin{Bmatrix} {{\overset{.}{\Gamma}}_{T}(t)} \\ {{\overset{.}{\Gamma}}_{W}(t)} \end{Bmatrix}} + {\begin{bmatrix} K_{\Gamma;T} & \text{ } \\ \text{ } & K_{\Gamma;W} \end{bmatrix}\begin{Bmatrix} {\Gamma_{T}(t)} \\ {\Gamma_{W}(t)} \end{Bmatrix}}} =} \\ {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{\begin{bmatrix} P_{T}^{T} \\ {- P_{W}^{T}} \end{bmatrix}{K_{c}\left( {t,j,k} \right)}\begin{matrix} \left\lbrack P_{T} \right. & {\left. {- P_{W}} \right\rbrack\begin{Bmatrix} {{\Gamma_{T}(t)} - {\Gamma_{T}\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \\ {{\Gamma_{W}(t)} - {\Gamma_{W}\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \end{Bmatrix}} \end{matrix}} + {\begin{bmatrix} P_{T}^{T} \\ {- P_{W}^{T}} \end{bmatrix}{F_{0}\left( {t,j,k} \right)}}} \right\}}} \end{matrix}$ where M_(Γ;T), C_(Γ;T) and K_(Γ;T) are modal mass, modal damping and modal stiffness matrices of cutting tool; M_(Γ;W), C_(Γ;W) and K_(Γ;W) are modal mass; modal damping and modal stiffness matrices of workpiece; {umlaut over (Γ)}_(T)(t), {dot over (Γ)}_(T) (t) and Γ_(T)(t) are modal acceleration, modal velocity and modal displacement vectors of cutting tool; {umlaut over (Γ)}_(W)(t), Γ_(W)(t) and Γ_(W)(t) are modal acceleration, modal velocity and modal displacement vectors of workpiece; P_(T), and P_(W) are mode shape matrices of cutting tool and workpiece, respectively; P_(T) ^(T) and P_(W) ^(T) are transposed matrices P_(T) and P_(W), respectively; for milling system with flexible cutting tool and rigid workpiece, the dynamic model becomes: $\begin{matrix} {{{M_{\Gamma;T}{\overset{¨}{\Gamma}}_{T}(t)} + {C_{\Gamma;T}{\overset{.}{\Gamma}}_{T}(t)} + {K_{\Gamma;T}\Gamma_{T}(t)}} =} \\ {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{P_{T}^{T}K_{c}\left( {t,j,k} \right)\left\{ {P_{T}\left\lbrack {{\Gamma_{T}(t)} - {\Gamma_{T}\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack} \right\}} + {P_{T}^{T}F_{0}\left( {t,j,k} \right)}} \right\}}} \end{matrix}$ for milling system with rigid cutting tool and flexible workpiece, the dynamic model becomes: $\begin{matrix} {{{M_{\Gamma;W}{\overset{¨}{\Gamma}}_{W}(t)} + {C_{\Gamma;W}{\overset{.}{\Gamma}}_{W}(t)} + {K_{\Gamma;W}\Gamma_{W}(t)}} =} \\ {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}{\left\{ {{P_{W}^{T}{K_{c}\left( {t,j,k} \right)}\left\{ {P_{W}\left\lbrack {{\Gamma_{W}(t)} - {\Gamma_{W}\left( {t - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack} \right\}} - {P_{W}^{T}{F_{0}\left( {t,j,k} \right)}}} \right\}.}}} \end{matrix}$
 3. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 2 comprises sub-steps of: step 2.1: performing state space transformation on the dynamic model in modal space by introducing state variable vector x(t)=[(Γ(t))^(T),({dot over (Γ)}(t))^(T)]^(T), and delay differential equation with multiple delays in state space becomes: ${\overset{\cdot}{x}(t)} = {{{Ax}(t)} + {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{{B\left( {t,j,k} \right)}\left\lbrack {{x(t)} - {x\left( {t - {\sum_{p = {\overset{\_}{k}}_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack} + {D\left( {t,j,k} \right)}} \right\}}}}$ where x(t) denotes state variable vector, {dot over (x)}(t) denotes first derivative of x(t) with regard to time t, ${A = \begin{bmatrix} 0 & I \\ {{- M_{\Gamma}^{- 1}}K_{\Gamma}} & {{- M_{\Gamma}^{- 1}}C_{\Gamma}} \end{bmatrix}},$ ${{B\left( {t,j,k} \right)} = \begin{bmatrix} 0 & 0 \\ {M_{\Gamma}^{- 1}P^{T}{K_{c}\left( {t,j,k} \right)}P} & 0 \end{bmatrix}},$ ${{D\left( {t,j,k} \right)} = \begin{bmatrix} 0 \\ {M_{\Gamma}^{- 1}P^{T}{F_{0}\left( {t,j,k} \right)}} \end{bmatrix}},$ and I is unit matrix; step 2.2: dividing one rotation period of cutting tool T into 2i_(m)+2 equal intervals with a series of discrete nodes {t₀,t₁, . . . t_(2i) _(m) ₊₂}; analytical expression of dynamic responses at interval [t_(i),t] with regard to the delay differential equation with multiple delays in step 2.1 is: ${x(t)} = {{e^{A({\tau - t_{i}})}{x\left( t_{i} \right)}} + {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}{\int_{t_{i}}^{t}{\left\{ {{e^{A({t - \xi})}{{B\left( {\xi,j,k} \right)}\left\lbrack {{x(\xi)} - {x\left( {\xi - {\sum_{p = k_{v}}^{k}{\tau\left( {j,p} \right)}}} \right)}} \right\rbrack}} + {e^{A({t - \xi})}{D\left( {\xi,j,k} \right)}}} \right\} d\xi}}}}}$ where ξ is integration variable; for brevity, define x_(i):=x(t_(i)), E(t,t_(i)):=e^(A(t−t) ^(i) ⁾x(t_(i)), H_(j,k)(t,ξ,x(ξ)):=e^(A(t−ξ))B(ξ,j,k)[x(ξ)−x(ξ−τ(j,k,p))] and J_(j,k)(t,ξ):=e^(A(t+ξ))D(ξ,j,k); on sub-interval [t_(2i),t_(2i+2)], x(t) is approximated using Simpson's rule: $x_{{2i} + 2} = {{E\left( {t_{{2i} + 2},t_{2i}} \right)} + {\frac{h}{3}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{H_{j,k}\left( {t_{{2i} + 2},t_{2i},x_{2i}} \right)} + {4{H_{j,k}\left( {t_{{2i} + 2},t_{{2i} + 1},x_{{2i} + 1}} \right)}} + {H_{j,k}\text{ }\left( {t_{{2i} + 2},t_{{2i} + 2},x_{{2i} + 2}} \right)} + {J_{j,k}\left( {t_{{2i} + 2},t_{2i}} \right)} + {4{J_{j,k}\left( {t_{{2i} + 2},t_{{2i} + 1}} \right)}} + {J_{j,k}\left( {t_{{2i} + 2},t_{{2i} + 2}} \right)}} \right\}}}}}$ on sub-interval [t_(2i),t_(2i+1)], x(i) is approximated using fourth-order Runge-Kutta algorithm: $x_{{2i} + 1} = {{E\left( {t_{{2i} + 1},t_{2i}} \right)} + {\frac{h}{6}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left\{ {{H_{j,k}\left( {t_{{2i} + 1},t_{2i},x_{2i}} \right)} + {2{H_{j,k}\left( {t_{{2i} + 1},t_{{2i} + {1/2}},x_{{2i} + {1/2}}} \right)}} + {2{H_{j,k}\left( {t_{{2i} + 1},t_{{2i} + {1/2}},x_{{2i} + {1/2}}} \right)}} + {H_{j,k}\left( {t_{{2i} + 1},t_{{2i} + 1},x_{{2i} + 1}} \right)} + {J_{j,k}\left( {t_{{2i} + 1},t_{2i}} \right)} + {2{J_{j,k}\left( {t_{{2i} + 1},t_{{2i} + {1/2}}} \right)}} + {2{J_{j,k}\left( {t_{{2i} + 1},t_{{2i} + {1/2}}} \right)}} + {J_{j,k}\left( {t_{{2i} + 1},t_{{2i} + 1}} \right)}} \right\}}}}}$ where state variable vector at middle points t_(2i+1/2) is approximated with Lagrange formula: x _(2i+1/2)=⅜x _(2i)+¾x _(2i−1)−⅛x _(2i+2) re-organizing the analytical expression of dynamic responses on sub-interval [t_(2i),t_(2i+2)] into: ${{F_{2i}x_{2i}} + {F_{{2i} + 1}x_{{2i} + 1}} + {F_{{2i} + 2}x_{{2i} + 2}}} = {{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left( {{F_{{2i} - {\tau({j,k})}}x_{{2i} - {\tau({j,k})}}} + {F_{{2i} + 1 - {\tau({j,k})}}x_{{2i} + 1 - {\tau({j,k})}}} + {F_{{2i} + 2 - {\tau({j,k})}}x_{{2i} + 2 - {\tau({j,k})}}}} \right)}} + {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left( {S_{2i} + S_{{2i} + {1/2}} + S_{{2i} + 1}} \right)}}}$ where ${F_{2i} = {{- e^{Ah}} - {\frac{h}{6}e^{Ah}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{2i},j,k}}}} - {{\frac{3}{8} \cdot \frac{4h}{6}}e^{{Ah}/2}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + {1/2}},j,k}}}}}},$ ${F_{{2i} + 1} = {I - {{\frac{3}{4} \cdot \frac{4h}{6}}e^{{Ah}/2}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + {1/2}},j,k}}}} - {\frac{h}{6}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + 1},j,k}}}}}},$ ${F_{{2i} + 2} = {{\frac{1}{8} \cdot \frac{4h}{6}}e^{{Ah}/2}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + {1/2}},j,k}}}}},$ ${F_{{2i} - {\tau({j,k,p})}} = {{{- \frac{h}{6}}e^{Ah}B_{{2i},j,k}} - {{\frac{3}{8} \cdot \frac{4h}{6}}e^{{Ah}/2}B_{{{2i} + {1/2}},j,k}}}},$ ${F_{{2i} + 1 - {\tau({j,k,p})}} = {{{{- \frac{3}{4}} \cdot \frac{4h}{6}}e^{{Ah}/2}B_{{{2i} + {1/2}},j,k}} - {\frac{h}{6}B_{{{2i} + 1},j,k}}}},$ ${F_{{2i} + 2 - {\tau({j,k,p})}} = {{\frac{1}{8} \cdot \frac{4h}{6}}e^{{Ah}/2}B_{{{2i} + {1/2}},j,k}}},$ ${S_{2i} = {\frac{h}{6}e^{A \cdot h}D_{{2i},j,k}}},$ ${S_{{2i} + {1/2}} = {{\frac{h}{6} \cdot 4}e^{{A \cdot h}/2}D_{{{2i} + {1/2}},j,k}}},$ ${S_{{2i} + 1} = {\frac{h}{6}D_{{{2i} + 1},j,k}}},$ and h is the discretization step; re-organizing the analytical expression of dynamic responses on sub-interval [t_(2i),t_(2i+2)] into: ${{G_{2i}x_{2i}} + {G_{{2i} + 1}x_{{2i} + 1}} + {G_{{2i} + 2}x_{{2i} + 2}}} = {{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left( {{G_{{2i} - {\tau({j,k})}}x_{{2i} - {\tau({j,k})}}} + {G_{{2i} + 1 - {\tau({j,k})}}x_{{2i} + 1 - {\tau({j,k})}}} + {G_{{2i} + 2 - {\tau({j,k})}}x_{{2i} + 2 - {\tau({j,k})}}}} \right)}} + {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\left( {T_{2i} + T_{{2i} + 1} + T_{{2i} + 2}} \right)}}}$ where ${G_{2i} = {{- e^{2{Ah}}} - {\frac{h}{3}e^{2{Ah}}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{2i},j,k}}}}}},$ ${G_{{2i} + 1} = {{- \frac{4h}{3}}e^{Ah}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + 1},j,k}}}}},$ ${G_{{2i} + 2} = {I - {\frac{h}{3}{\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}B_{{{2i} + 2},j,k}}}}}},$ ${G_{{2i} - {\tau({j,k,p})}} = {{- \frac{h}{3}}e^{2Ah}B_{{2i},j,k}}},$ ${G_{{2i} + 1 - {\tau({j,k,p})}} = {{- \frac{4h}{3}}e^{Ah}B_{{{2i} + 1},j,k}}},$ ${G_{{2i} + 2 - {\tau({j,k,p})}} = {{- \frac{h}{3}}B_{{{2i} + 2},j,k}}},$ ${T_{2i} = {\frac{h}{3}e^{{A \cdot 2}h}D_{{2i},j,k}}},$ $T_{{2i} + 1} = {{\frac{h}{3} \cdot 4}e^{A \cdot h}D_{{{2i} + 1},j,k}}$ and ${T_{{2i} + 2} = {\frac{h}{3}D_{{{2i} + 2},j,k}}};$ step 2.3: based on the formulae in step 2.2, establishing discrete mapping relationship between two adjacent rotation periods of cutting tool [−T,0] and [0,T]: P₁y_([0, T]) = Qy_([−T, 0]) + P₂y_([0, T]) + z_([0, T]) where ${y_{\lbrack{0,T}\rbrack} = \begin{bmatrix} x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\  \vdots \\ x_{{2i_{m}} + 1} \\ x_{{2i_{m}} + 2} \end{bmatrix}},{y_{\lbrack{{- T},0}\rbrack} = \begin{bmatrix} x_{0 - T} \\ x_{1 - T} \\ x_{2 - T} \\ x_{3 - T} \\ x_{4 - T} \\  \vdots \\ x_{{2i_{m}} + 1 - T} \\ x_{{2i_{m}} + 2 - T} \end{bmatrix}},{z_{\lbrack{0,T}\rbrack} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\begin{bmatrix} 0 \\ {S_{0} + S_{1/2} + S_{1}} \\ {T_{0} + T_{1} + T_{2}} \\ {S_{2} + S_{2 + {1/2}} + S_{3}} \\ {T_{2} + T_{3} + T_{4}} \\  \vdots \\ {S_{2i_{m}} + S_{{2i_{m}} + {1/2}} + S_{{2i_{m}} + 1}} \\ {T_{2i_{m}} + T_{{2i_{m}} + 1} + T_{{2i_{m}} + 2}} \end{bmatrix}}}},$ ${P_{1} = \begin{bmatrix} I & 0 & 0 & 0 & 0 & \ldots & 0 & 0 & 0 \\ F_{0} & F_{1} & F_{2} & 0 & 0 & \ldots & 0 & 0 & 0 \\ G_{0} & G_{1} & G_{2} & 0 & 0 & \ldots & 0 & 0 & 0 \\ 0 & 0 & F_{2} & F_{3} & F_{4} & \ldots & 0 & 0 & 0 \\ 0 & 0 & G_{2} & G_{3} & G_{4} & \ldots & 0 & 0 & 0 \\  \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & 0 & \ldots & F_{2i_{m}} & F_{{2i_{m}} + 1} & F_{{2i_{m}} + 2} \\ 0 & 0 & 0 & 0 & 0 & \ldots & G_{2i_{m}} & G_{{2i_{m}} + 1} & G_{{2i_{m}} + 2} \end{bmatrix}},$ $Q = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}\begin{bmatrix} 0 & \ldots & 0 & 0 & 0 & \ldots & 0 & 0 & {I/{NN}_{a}} \\ 0 & \ldots & F_{- {\tau({j,k,p})}} & F_{1 - {\tau({j,k,p})}} & F_{2 - {\tau({j,k,p})}} & \ldots & 0 & 0 & 0 \\ 0 & \ldots & G_{- {\tau({j,k,p})}} & G_{1 - {\tau({j,k,p})}} & G_{2 - {\tau({j,k,p})}} & \ddots & \vdots & 0 & 0 \\  \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & \ldots & 0 & 0 & 0 & \ldots & F_{{2i} - {\tau({j,k,p})}} & F_{{2i} + 1 - {\tau({j,k,p})}} & F_{{2i} + 2 - {\tau({j,k,p})}} \\ 0 & \ldots & 0 & 0 & 0 & \ldots & G_{{2i} - {\tau({j,k,p})}} & G_{{2i} + 2 - {\tau({j,k,p})}} & G_{{2i} + 2 - {\tau({j,k,p})}} \\  \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & \ldots & 0 & 0 & 0 & \ldots & 0 & 0 & 0 \end{bmatrix}}}$ and $P_{2} = {\sum\limits_{k = 1}^{N}{\sum\limits_{j = 1}^{N_{a}}{\begin{bmatrix} 0 & 0 & 0 & \ldots & 0 & 0 & 0 & \ldots & 0 \\  \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ldots & 0 \\ F_{{2i} + 2 + {\tau({j,k,p})}} & F_{{2i} + 3 + {\tau({j,k,p})}} & F_{{2i} + 4 - {\tau({j,k,p})}} & \ldots & 0 & 0 & 0 & \ldots & 0 \\ G_{{2i} + 2 - {\tau({j,k,p})}} & G_{{2i} + 3 - {\tau({j,k,p})}} & G_{{2i} + 4 - {\tau({j,k,p})}} & \ddots & 0 & 0 & 0 & \ldots & 0 \\ 0 & 0 & \vdots & \ddots & 0 & 0 & 0 & \ldots & 0 \\  \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \ldots & F_{{2i_{m}} - {\tau({j,k,p})}} & F_{{2i_{m}} + 1 - {\tau({j,k,p})}} & F_{{2i_{m}} + 2 - {\tau({j,k,p})}} & \ldots & 0 \\ 0 & 0 & 0 & \ldots & G_{{2i_{m}} - {\tau({j,k,p})}} & G_{{2i_{m}} + 1 - {\tau({j,k,p})}} & G_{{2i_{m}} + 2 - {\tau({j,k,p})}} & \ldots & 0 \end{bmatrix}.}}}$
 4. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 3 is: constructing transition matrix Φ of the milling system: $\Phi = \left\{ \begin{matrix} {{\left( {P_{1} - P_{2}} \right)^{- 1}Q},} & {{{if}{❘{P_{1} - P_{2}}❘}} \neq 0} \\ {{\left( {P_{1} - P_{2}} \right)^{\dagger}Q},} & {{{if}{❘{P_{1} - P_{2}}❘}} = 0} \end{matrix} \right.$ where |P₁−P₂| denotes determinant of matrix (P₁−P₂); † denotes Moore-Penrose pseudoinverse of matrix (P₁−P₂); based on Floquet theory, stability of milling system depends on eigenvalue of Φ:it is stable if all modules of eigenvalues of Φ are less than 1; it is unstable otherwise.
 5. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 4 is: based on fixed mapping point theorem, the state variable vector x(t) satisfies x(t_(i))=x(t_(i)−T) for chatter-free cutting conditions; therefore, it has y_([0,T])=y_(−T,0); calculating state variables on all discrete points by: $y^{*} = {y_{\lbrack{0,T}\rbrack} = {y_{\lbrack{{- T},0}\rbrack} = \left\{ \begin{matrix} {{\left( {P_{1} - P_{2} - Q} \right)^{- 1}z_{\lbrack{0,T}\rbrack}},} & {{{if}{❘{P_{1} - P_{2} - Q}❘}} \neq 0} \\ {{\left( {P_{1} - P_{2} - Q} \right)^{\dagger}z_{\lbrack{0,T}\rbrack}},} & {{{if}{❘{P_{1} - P_{2} - Q}❘}} = 0} \end{matrix} \right.}}$ where y* consists of modal displacement vector Γ(t_(i)) and modal velocity vector {dot over (Γ)}(t_(i)) at all discrete time nodes i=0,1, . . . ,2i_(m)+2 during one rotation period of cutting tool.
 6. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 5 comprises sub-steps of: step 5.1: transforming modal displacements into physical space, and calculating relative vibration displacements q(t_(i)) between cutting tool and workpiece: q(t _(i))=q _(T)(t _(i))−q _(W)(t _(i))=P _(T)Γ_(T)(t _(i))−P _(W)Γ_(W)(t _(i)) where q_(T) (t_(i)) and q_(W)(t_(i)) are vibration displacements of cutting tool and workpiece, respectively; step 5.2: extracting vibration displacements in direction normal to milled surface of workpiece (defined as {right arrow over (Oy)} direction) from q(t_(i)), which constructs a new vector Q*_(y): Q _(y)*={[q _(y)(t ₁)]T,[q _(y)(t ₂)]^(T), . . . ,[q _(x)(t _(2i) _(m) ₊₂)]T} extracting vibration displacements in feed direction (defined as {right arrow over (Ox)} direction), which constructs a new vector Q_(y)*: Q _(x)*={[q _(x)(t ₁)]^(T),[q _(x)(t ₂)]^(T), . . . ,[q _(x)(t _(2i) _(m) ₊₂)]T} step 53: based on kinematic synthesis of relative feed movement between cutting tool and workpiece, rotation of cutting tool and relative vibrations between cutting tool and workpiece, calculating motion trajectories in {right arrow over (Oy)} direction and {right arrow over (Ox)} direction for cutting element (j,k), respectively: ${\Theta_{y}^{*}\left( {i,j,k} \right)} = {\underset{feed}{\underset{︸}{\frac{f \cdot t_{i}}{1000 \times 60}{\sin\left( \varphi_{fy} \right)}}} - \underset{rotation}{\underset{︸}{{R\left( {j,k} \right)}{\cos\left( {\varphi_{ref} - {\phi\left( {j,k} \right)} - \frac{z_{j}\tan\left( \beta_{k} \right)}{R} + {2\pi\frac{\Omega t_{i}}{60}}} \right)}}} + \underset{vibration}{\underset{︸}{Q_{y}^{*}\left( {i,j} \right)}}}$ ${\Theta_{y}^{*}\left( {j,k,i} \right)} = {\underset{feed}{\underset{︸}{\frac{f \cdot t_{i}}{1000 \times 60}{\cos\left( \varphi_{fy} \right)}}} + \underset{rotation}{\underset{︸}{{R\left( {j,k} \right)}\sin\left( {\varphi_{ref} + {\phi\left( {j,k} \right)} - \frac{z_{j}\tan\left( \beta_{k} \right)}{R} + \frac{60t_{i}}{\Omega}} \right)}} + \underset{vibrations}{\underset{︸}{Q_{x}^{*}\left( {j,i} \right)}}}$ where f is relative feed rate between cutting tool and workpiece in mm/min; φ_(fy) is angle between feed direction of cutting tool and normal direction of workpiece; R(j,k) is actual cutting radius of cutting element (j,k); φ_(ref) is reference angle used in the Generalized Runge-Kutta method; ϕ(j,k) is pitch angle between Tooth k and Tooth 1 for j-th discrete axial slice; β_(k) is helix angle with regard to Tooth k; z_(f) is axial height of j-th discrete axial slice; R is geometric radius of cutting tool; Ω is spindle rotation speed; Q_(y)*(j,i) and Q_(x)*(j,i) are vibration displacements of cutting element (j,k) in {right arrow over (Oy)} and {right arrow over (Ox)} directions, respectively.
 7. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 6 comprises sub-steps of: step 6.1: selecting part of motion trajectories of cutting tool that are near milled surface of workpiece based on following equation, which construct motion trajectory points to be densified with interpolation, i.e., (Θ_(x_trim)*(j,k,i), Θ_(y_trim)*(j,k,i)): Θ_(y)*(j,k,i)≥max{Θ_(y)*(j,k,i)}−αR, down-milling Θ_(y)*(j,k,i)≤min{(Θ_(y)*(j,k,i)}+αR, up-milling where α is a coefficient to adjust selecting range; step 6.2: determining lower and upper limits of the selected trajectory points in {right arrow over (Ox)} direction, i.e., x_(min)=min{Θ_(x_trim)*(j,k,i)} and x_(max)=max{Θ_(x_trim)*(j,k,i)}; discretize interval [x_(min),x_(max)] equally with a step δx to obtain x coordinate set of interpolation points {x_(min), x_(min)+δx,x_(min)+2·δx, . . . ,x_(max)}, where number of points is denoted by N_(s); step 6.3: taking the selected points (Θ_(x_trim)*(j,k,i), Θ_(y_trim)*(j,k,i)) in step 6.1 as known points, using spline interpolation algorithm to determine y coordinates of interpolation points with regard to x coordinates {x_(min), x_(min)+δx,x_(min)+2·δx, . . . ,x_(max)}, which results in densified motion trajectories of cutting tool edges during one rotation period T of cutting tool, i.e., (x_(s)(l), y_(s)(l)),l=1, 2 . . . N_(s); step 6.4: based on periodicity of chatter-free milling dynamic responses, i.e., x_(s)(l) and x_(s)(l)+n_(rev)T share same y coordinate y_(s)(l), obtaining densified motion trajectory points of cutting tool during n_(rev) rotation periods of cutting tool, i.e., (x_(s_n)(l), y_(s_n)(l)), l=1, 2, . . . N_(s_n), where N_(s_n), is number of interpolation points on n_(rev) rotation periods of cutting tool, and n_(rev)≥4; step 6.5: selecting densified motion trajectory points of cutting tool edges on middle periods that satisfy x_(min)+T≤x_(s_n)(l)≤n_(rev)·x_(max)−T to construct a new densified point set (x_(s_n_trim)(l), y_(s_n_trim)(l), l=1, 2, . . . N_(s_n_trim), where N_(s_n_trim) is number of trimmed interpolation points; step 6.6: sorting the selected set (x_(s_n_trim)(l), y_(s_n_trim)(l)), in order of axial height to construct N_(a) sub-sets of densified points, i.e., (x_(s_n_trim)(j,l), y_(s_n_trim)(j,l)), j=1,2, . . . ,N_(a); step 6.7: at each axial height, comparing values of v coordinates for all teeth that have same x coordinate; based on following equation, selecting densified points that are nearest to milled surface of workpiece to construct final surface topography points of workpiece (x_(surf)(j,l), y_(surf) (j,l)), l=1, 2, . . . N_(s_n_trim): $\begin{matrix} {{{y_{surf}\left( {j,l} \right)} = {\max\limits_{k}\left\{ {y_{{s\_ n}{\_ trim}}\left( {j,l} \right)} \right\}}},} & {{down} - {milling}} \\ {{{y_{surf}\left( {j,l} \right)} = {\min\limits_{k}\left\{ {y_{{s\_ n}{\_ trim}}\left( {j,l} \right)} \right\}}},} & {{up} - {milling}} \end{matrix}$
 8. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 7 is: calculating surface location error (SLE) based on motion trajectories Θ_(y)*(i,j,k) in normal direction of workpiece; surface location error SLE(j,k) with regard to cutting element (j,k) is: ${{SLE}\left( {j,k} \right)} = \left\{ \begin{matrix} {{{\max\limits_{i}\left\{ {{\Theta_{y}^{*}\left( {j,k,i} \right)}{❘{1 \leq i \leq {{2i_{m}} + 2}}}} \right\}} - {D/2}},} & {{down} - {milling}} \\ {{{\underset{i}{- \min}\left\{ {\Theta_{y}^{*}\left( {j,k,i} \right){❘{1 \leq i \leq {{2i_{m}} + 2}}}} \right\}} - {D/2}},} & {{up} - {milling}} \end{matrix} \right.$ surface location error with regard to axial height z_(j) is: ${{SLE}(j)} = \left\{ \begin{matrix} {{{\max\limits_{i,k}\left\{ {{{\Theta_{y}^{*}\left( {j,k,i} \right)}❘{1 \leq i \leq {{2i_{m}} + 2}}},{1 \leq k \leq N}} \right\}} - {D/2}},} & {{down} - {milling}} \\ {{{\underset{i,k}{- \min}\left\{ {{{\Theta_{y}^{*}\left( {j,k,i} \right)}❘{1 \leq i \leq {{2i_{m}} + 2}}},{1 \leq k \leq N}} \right\}} - {D/2}},} & {{up} - {milling}} \end{matrix} \right.$ where positive value mean overcut, and negative value mean undercut.
 9. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 8 is: based on points (x_(surf)(j,l), y_(surf)(j,l), y_(surf)(j,l)) that participate in generating surface topography of workpiece, calculating surface roughness: ${{Ra}(j)} = {\frac{1}{N_{{s\_ n}{\_ trim}}}{\sum\limits_{i^{\prime} = 1}^{N_{{s\_ n}{\_ trim}}}{{❘{{y_{surf}\left( {j,l} \right)} - {{\overset{\_}{y}}_{surf}\left( {j,l} \right)}}❘}.}}}$ 